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C. ORTNER 

Abstract. For a general class of atomistic-to-contirmum coupling methods, coupling multi- 
body interatomic potentials with a Pi-finite element discretisation of Cauchy-Born nonlinear 
elasticity, this paper adresses the question whether patch test consistency (or, absence of ghost 
forces) implies a first-order error estimate. 

In two dimensions it is shown that this is indeed true under the following additional technical 
assumptions: (i) an energy consistency condition, (ii) locality of the interface correction, (iii) 
volumetric scaling of the interface correction, and (iv) connectedness of the atomistic region. 
The extent to which these assumptions are necessary is discussed in detail. 



1. Introduction 

Defects in crystalline materials interact through their elastic fields far beyond their atomic 
neighbourhoods. An accurate computation of such defects requires the use of atomistic mod- 
els; however, the size of the atomistic systems that are often required to accurately represent 
the elastic far-field makes atomistic models infeasible or, at the very least, grossly inefficient. 
Indeed, atomistic accuracy is often only required in a small neighbourhood of the defect, while 
the elastic far field may be approximated using an appropriate continuum elasticity model. 

Atomistic-to-continuum coupling methods (a/c methods) aim to exploit this fact by retaining 
atomistic models only in small neighbourhoods of defects, and coupling these neighbourhoods to 
finite element discretisations of continuum elasticity models; see Figure []Jc,d). By employing 
a coarse discretisation of the continuum model, such a process can achieve a considerable 
reduction in computational complexity, however, some of the first a/c methods suffered from 
the so-called "ghost force problem" : While homogeneous deformations are equilibria of both 
the pure atomistic and the pure continuum model, they are not equilibria of certain a/c models 
|37[ 150] [35] due to spurious forces — the "ghost forces" — that can arise at the interface between 
the atomistic and continuum regions. 

Much of the literature on a/c methods has focused on constructing a/c methods that did not 
exhibit, or reduced the effect of the "ghost forces" Ell E21 [13 SI Effl [13 E31 E2] • 

A straightforward solution was the introduction of force-based (i.e, non-conservative) methods 
[25"1 |5"U1 |5"T] [TP] [2U1 [33]. The construction of accurate energy-based coupling mechanisms turned 
out to be more challenging. Several creative approaches providing partial solutions to the 
problem were suggested [52l [151 SHI 122] , however, no general solution exists so far. 

The inconsistency of early a/c methods is reminiscent of the inconsistency problems encoun- 
tered in the early history of finite element methods. A simple criterion to test consistency of 
finite element methods is the patch test introduced by Irons et al [5]; see also [531 E]- The 
"ghost force problem" discussed above is precisely the failure of such a patch test. It is well 
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known that, in general, the patch test is neither necessary nor sufficient for convergence of 
finite element methods; see, e.g., [531 [6] where several variants of patch tests are discussed. The 
same is true for a/c methods: It was shown in |31j that a particular flavour of force-based a/c 
coupling typically has a consistency error of nearly 100%, even though it does pass the patch 
test. 

Although a growing numerical analysis literature exists on the subject of a/c methods |28 } I42 | 
[TOt [T21 HUl [T3l 1361 [TBI 130] , it has so far focused primarily on one- dimensional model problems. 
(A notable exception is the work of Lu and Ming on force-based hybrid methods [30]. However, 
the techniques used therein require large overlaps and cannot accommodate sharp interfaces.) 
Only specific methods are analyzed; the question whether absence of "ghost forces" (or, patch 
test consistency, as we shall call it) in general implies satisfactory accuracy has neither been 
posed nor addressed so far. The purpose of the present paper is to fill precisely this gap. 
After introducing a general atomistic model and a general class of abstract a/c methods, and 



establishing the necessary analytical framework, it will be shown in Theorem 10, which is the 
main result of the paper, that in two dimensions patch test consistency together with additional 
technical assumptions implies first-order consistency of energy-based a/c methods. 

1.1. Outline. jj2] gives a detailed introduction to the construction of a/c methods. This section 
also develops a new notation that is well suited for the analysis of 2D and 3D a/c methods. In 



£2.3.3 we give a precise statement of the patch test consistency condition. 

^3j contains a general framework for the a priori error analysis of a/c methods in W 1 '^- 
norms, similar to an error analysis of Galerkin methods with variational crimes. Several new 
technical results are presented in this section, such as the introduction of an oscillation operator 
to measure local smoothness of discrete functions ( ^3.2.1 ), an interpolation error estimate for 



piecewise affine functions ( ^3.2.2 ), and making precise the assumption made in much of the 



a/c numerical analysis literature that it is sufficient to estimate the modelling error without 
coarsening (£ |3.5[ ), The main ingredient left open in this analysis is a stability assumption, 
which requires a significant amount of additional work and is beyond the scope of this paper. 

§|4] presents two ID examples for modelling error estimates, which motivate the importance of 
patch test consistency, and discusses the modelling error estimates that can at best be expected 
if a method is patch test consistent. 

^5] introduces the two main auxiliary results used in the 2D analysis of (i) Shapeev's bond 
density lemma, which allows the translation between bond integrals and volume integrals; and 
(ii) a representation theorem for discrete divergence-free Po-tensor fields. 

Finally, ^6] establishes the main result of this paper, Theorem 10 if an a/c method is 



patch test consistent and satisfies other natural technical assumptions, then it is also first-order 
consistent. The proof depends on a novel construction of stress tensors for atomistic models, 
related to the virial stress [2] (generalizing the ID construction in [32^133]). and a corresponding 
construction for the stress tensor associated with the a/c energy. Moreover, we discuss in detail 



to what extent the technical conditions of Theorem 10 are required. For example, we show 
that, if the atomistic region is finite (i.e., completely surrounded by the continuum region) then 
the condition of global energy consistency already follows from patch test consistency. 

1.2. Sketch of the main result. In this section we discuss the main result in non-rigorous 
terms. Let be an atomistic energy functional and let <#ac be an a/c energy functional, which 
uses the atomistic description in part of the computational domain, and couples it to a finite 
element discretisation of Cauchy-Born nonlinear elasticity (cf. ^2.3.2 for its definition), with 



suitable interface treatment. Suppose that the following conditions are satisfied: 

(i) <f ac is patch test consistent: Every homogeneous deformation is a critical point of 
(cf- §0- 
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(ii) is globally energy consistent: <§& c is exact for homogeneous deformations (cf. § 2.3.2). 

(iii) Locality and scaling of the interface correction: The interface correction has (roughly) 
the same interaction range as the atomistic model, and has volumetric scaling (cf. 



(iv) The atomistic region f2 a is connected. 

(v) Stability: For some p G [l>oo], and for deformations y in a neighbourhood of the 
atomistic solution, the second variation <5 2 (=? ac (y) is stabl e, w hen understood as a linear 



operator from (discrete variants of) W 1,p to W ' p (cf. ^3.3). 
(vi) The finite element mesh in the continuum region is shape regular. [7J 

Under conditions (v) and (vi), we show in ^3] that the error between the atomistic solution 
y a and the a/c solution y ac can be bounded by 

||Vy a - Vy ac || LP(Q) < £r dcl + £t + HhVVII^, (1) 

where Q is the computational domain, Q c the continuum region, h a local mesh size function, 
£^ xt is the consistency error for the treatment of external forces, and £;™ odel is the modelling 
error, which we describe next. Since y a is not a continuous deformation, but a deformation of a 
discrete lattice, the various terms appearing above need to be interpreted with care. This is the 
focus of ^3j For example, in the rigorous version of ([!]), we will replace V 2 y a by an oscillation 
of a suitably defined gradient. 

The step outlined above does not distinguish different variants of a/c methods. The error 
introduced by the coupling mechanism and the continuum model is contained in the modelling 
error 

gmodel = ||^ ac(ya) _ ^ a (y a )|| w _ lip) 

where W<T 1,P is a suitably defined negative Sobolev norm on an atomistic grid. If <f ac is not 
patch test consistent, then, typically, £™ odel = O(l). However, if conditions (i)-(iv) hold, and 
if the problem is set in either one or two space dimensions, then we will prove in ^6] that 

pmodel < ||V72 II ^ 

where e is the atomistic spacing and f^i an interface region. A rigorous statement of ([2]) is given 



in Theorem 10 which is the main result of the paper. 



1.3. Basic notational conventions. Vectors are denoted by lower case roman symbols, x,y, 
a, b, and so forth. Matrices are denoted by capital sans serif symbols A, B, F, G, and so forth. 
We will not distinguish between row and column vectors. If two vectors are multiplied, then it 
will be specified whether the operation is the dot product or the tensor product: for a, b 6 R fc , 
we define 

k 

a • b := 2. a jbj, and a<g> b := (aj6j)i = i k, 

where, in a <g> b, i is the row index and j the column index. 

When matrix fields take the role of stress tensors, we will also call them tensor fields and use 
greek letters a, E, and so forth. 

The Euclidean norm of a vector and the Frobenius norm of a matrix are denoted by | • | . The 
f-norms, p £ [l,oo], of a vector or matrix are denoted by | • | p , and sometimes by || • For 
e > 0, the weighted f p -norms, on an index set oS^ 7 , are defined as 

ll a lk(j^) := £l/p \\ a \\ep(,y)- 
The topological dual of a vector space % is denoted by with duality pairing (•,-). 
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If A C M. d is a measurable set then \A\ denotes its d-dimensional volume. The symbol f A 
denotes J A , provided that \A\ > 0. If A has a well-defined area or length, then these are 

denoted, respectively, by area(A) and length^). 

For a measurable function / : A — > R, ||/||lp(j4) denotes the standard L p -norm. If / : A — > 
R k * m , then \\f\\ LP{A) ~ || \f\ p \\ LP{A) . 

Partial derivatives with respect to a variable Xj, say, are denoted by d/dxj. The Jacobi 
matrix of a differentiate function / : A — > M fc is denoted by df. The symbol d will also be 
redefined in some contexts, but used with essentially the same meaning as here. 

When / is a deformation or a displacement, then we will also write V/ = df. For r G W*, 
the uni-directional derivative is denoted by 

g/ W = Urn + *>-■»*>, 

K ' t\0 t 

whenever this limit exists. The symbol D r denotes a finite difference operator, which will be 
defined in pX2] 



Additional notation will be defined throughout. A list of symbols, with references to their 
definitions is given in ^Bj 

2. Introduction to Atomistic/Continuum Model Coupling 

In this section we introduce a general multi-body interaction model with periodic boundary 
conditions. This choice of boundary condition is crucial since the non-locality of the atomistic 
interactions makes an analysis of Dirichlet or Neumann boundaries challenging. The periodic 
boundary conditions can also be understood as "artificial boundary conditions" for infinite 
crystals. 

Next, we describe the construction of energy-based a/c methods that couple the atomistic 
interaction potential with a Pi-finite element discretization of the Cauchy-Born continuum 
model. We motivate the patch test ("ghost forces"), and why an interface correction is required 
to obtain accurate coupling schemes. 

2.1. An atomistic model with periodic boundary condition. 

2.1.1. Periodic deformations. Let d G {1,2,3} denote the space dimension. We will, in sub- 
sequent sections, restrict our analysis to d G {1, 2}, however, the introduction to a/c coupling 
methods is independent of the dimension. 

For some N £N, and e := 1/N, we define the periodic reference cell 

jgf : =e{ - N + l,...,N} d . 

The space of 2Z d -periodic displacements of Jzf* := eTL d is given by 

<% := {u : if # -»■ R d | u(x + £) = u(x) for all f G 2Z d ,x G JSf}. 

A homogeneous deformation (or, Bravais lattice) of Jzf* is a map : -> R d , where 

A G M. dxd and y^{x) = IKx for all x G Jzf#. We denote the space of periodic deformations of 
by 

^ := {y : Jz^ # -> M d | y - y k G W for some A G R dxd }, 
and the space of deformations with prescribed macroscopic strain A by 

W is a linear space, while is an affine subspace of & '. 
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For future reference we define notation that extends sets periodically: For any set A C M d 
we define A& := U^e2Z d (^ + ^)- This notation is consistent with the definition of Jzf*. If is 
a family of sets, then := {A# \ A G 



2.1.2. The atomistic energy. For a map v : Jz?* — > G N, and r G Z d \ {0}, we define the 
finite difference operator 

£Kz) := V{x + £r) e ~ v{x \ 
We assume that the stored elastic energy of a deformation y G ^ is given in the form 

^(y)--=e d ^2v(D^)), (3) 

where # C Z d \ {0} is a finite interaction range, Dcgy(x) := (D r y(x)) r& ^, and where V G 
(^((M^) 1 ^) is a multi-body interaction potential. Under these conditions, <f a G C 2 (^). 

The scaling of the lattice, of the finite difference operator, and of the energy were chosen to 
highlight the natural connection between molecular mechanics and continuum mechanics. For 
example, e d ^2 resembles an integral (or, a Riemann sum), while D r resembles a directional 
derivative. It should be stressed, however, that e is a fixed parameter of the problem, which is 
small but does not tend to zero. 

The formulation ^ includes all commonly employed classical interatomic potentials (see, 
e.g., |19} I47j): pair potentials such as the Lennard- Jones or Morse potential, bond-angle po- 
tentials, embedded atom potentials, bond-order potentials, or any combination of the former, 
provided that they have a finite interaction range. The generality of the interaction potential 
also includes effective potentials obtained for in-plane or anti-plane deformations of 3D crystals. 

No major difficulties should be expected in generalizing the analysis to infinite interaction 
ranges, provided the interaction strength decays sufficiently fast. A generalization of the anal- 
ysis to genuine long-range interactions such as Coulomb interactions is not obvious. 

2.1.3. Assumptions on the interaction potential. For g = (g r ) r &M £ (M. d )^ we denote the first 
and second partial derivatives of V at g, respectively, by 

^(g):=Pel d , and d r>s V(g) := ER dxd , for r, s, G St. 

dg r dg r dg s 

Throughout this work we assume the following global bound: there exist constants M r a s > 0, 
r, s G such that 

sup ||d r>s y(g)|| < M r a s , for all r, s G Si, (4) 
where || • || denotes the £ 2 -operator norm. This assumption contradicts realistic interaction 



models and is made to simplify the notation; see £2.2.1 for further discussion of this issue. 

In the subsequent analysis we will, in fact, never make direct use of the second partial 
derivatives d r ^ s V directly, but only use the resulting Lipschitz property, 

\d r V(g) - d r V(h)\ < Ks \Ss - h.\ Vg, h G (M d )'«, r G SI. (5) 

The proof is straightforward. For future reference, we also define the constant 

M a := Y, VW s \ M t,s- (6) 
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2.1.4. The variational problem. Let G C 2 (^ ;M) be the potential of external forces mod- 
elling, for example, a substrate or an indenter. As explained in §2.2.3[ one may also use such 
a potential to model simple point defects such as vacancies or impurities. 

Given a potential of external forces <^ a and a macroscopic strain A, we consider the problem 
of finding local minimizers of the total energy S^ ot := <§ & + & & in in short, 



'a 

atot / 



y a G argmm<C(2/), (7) 

where argmin denotes the set of local minimizers. 

If y a solves ([7]), then y a is a critical point of ^ a ot , that is, 

(S<?M + 5^a(2/a), n> = Vu G <2T, (8) 

where, for a functional <<? G C 1 (^), we define the /irsi variation of J 5 at y, as 

(6£(y),z) := ^(y + tz) \ t=0 for y, z G 
If <^ G C 2 (^) then the second variation is defined analogously as 

(5 2 g(y)zi,z 2 ) := —(S^(y + tz x ), z 2 )\ t=0 for y,zi,z 2 G ^. 

The same notation will be used for functionals defined on different spaces. 

2.1.5. The patch test for the atomistic model. The following proposition can be understood 
as the patch test for the atomistic model: In the absence of external forces and defects a 
homogeneous lattice is always a critical point of # a . 

Proposition 1. (5<f a (y A ), u) = for all u £ W and for all A G R dxd . 
Proof. Let y, z G then 

(5&(y), z) =e d Y,Y. drV{D M y(x)) ■ D r z(x). (9) 

Fix A G R dxd and u G ^, then 

(*<&(va),«> = E £<i E ^(A^) • 



Y d r V(A^) -\e d Y^ D r u{x)\ = 0, 



where we have used the fact that u is periodic. Above, and throughout, we use the notation 
Mi = (Ar) re # = D m y k {x). □ 

2.2. Remarks on the atomistic model. 



2.2.1. Invertibility of deformations. In i 2.1.2 and £2.1.3 we have assumed that S a is twice dif- 
ferentiable at all deformations y G & , and that the second partial derivatives of the interaction 
potential are globally bounded. 

However, realistic interaction potentials V take the value +oo if two atoms occupy the same 
position in space, and hence can only be differentiable at deformations that are one-to-one 
(i.e., "true" deformations). With only minor additional technicalities such potentials can be 
admitted in the analysis. The global bound Q would then be replaced by a local bound and 
certain explicit bounds on D&y; see, e.g., [4"0l HI]. 
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2.2.2. Reference cutoff. Another aspect of the atomistic energy ([3]), which makes it inappro- 
priate for realistic applications is that the interaction potential V has a cut-off radius in the 
reference configuration. In atomistic models, atoms are unconstrained in their position and 
hence two atoms that are far apart in the reference configuration may be arbitrarily close, and 
hence interact, in the deformed configuration. 

The reference cutoff in ^ is assumed only for the sake of brevity of the notation. One 



can, similarly as discussed in S 2.2.1 take a more general form of the interaction potential 
that does not suffer from this drawback and make suitable assumptions on deformations under 
consideration that control the interaction neighbourhood. 

2.2.3. Modelling crystal defects. Some simple crystal defects can be modelled via the poten- 
tial of the external forces, The simplest example is an impurity, where a single atom is 
replaced with an atom from a different species. For ([3]) this means that the interaction po- 
tential is changed from V(D^y(x)) to V mod (x; D^y(x)) in a neighbourhood of the impurity. 
Alternatively, one may keep the original form of <£* a and define 

= e d Yl [V mod & DavW) ~ V(D M y(x))] . 

Similarly, a vacancy can be modelled by simply removing all interactions with a given atom. 
This would yield a difficulty with the "unused" degrees of freedom for the position of the vacancy 
atom, which could simply be removed from the system [H]. An interstitial (an additional atom) 
can also be modelled fairly easily, but one would need to augment the variable y by additional 
degrees of freedom for the position of the interstitial atom. 

Dislocations, which possibly represent the most important class of crystal defects are, in 
general, more difficult to describe. In the atomistic minimization problem ([7]) they simply 
represent a special class of local minimizers, however, in the coupled atomistic/continuum 
models we discuss below most classes of dislocations less straightforward to embed; however, 
see [Ml El] for simple examples. 

2.3. Construction of a/c coupling methods. The atomistic model problem ([7]) is a finite- 
dimensional optimisation problem and is therefore, in principle, solvable using standard op- 
timisation algorithms. However, typical applications where atomistic models are employed 
require of the order 10 9 to 10 12 atoms or more |35} I37j. It is therefore desirable to construct 
computationally efficient coarse grained models. 

2.3.1. Galerkin projection. To motivate the idea of atomistic-to-continuum coupling we consider 
a crystal with a localized defect. Figure [l|a) shows a deformed 2D crystal with an impurity 
that repels its neighbouring atoms, causing a large local deformation. We observe that, except 
in a small neighbourhood of the defect, the atoms are arranged as a "smooth" deformation of 
the reference lattice Jzf*. It is therefore possible to approximate the atomistic configurations 
from a low-dimensional subspace constructed, for example, using a Pi-finite element method. 

Let CI := (— 1, and let ^ be a regular [7] triangulation of with vertices belonging 
to J2? # , that can be extended periodically to a regular triangulation Sffi of W 1 . We make the 
convention that elements T G 2?h are closed sets. For T G 3?h we define Ht '■= diam(T), and 
for each x G M. d we define h(x) := maxj/ij- 1 T G 2?h s.t. x G T}. 

We define the Pi finite element space 

Pi(^f) := {v h : -)• R | v h is piecewise affine w.r.t. ^ # }, 
and we denote the spaces of piecewise affine displacements and deformations, respectively, by 
^:=^nPi(5f ) d , %:=^nP!(5f ) d , and % := f A nPi(^ # ) d . 
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(C) 



(d) 



Figure 1. (a) A 2D atomistic configuration with an impurity that causes a 
large local deformation from the reference lattice el? . (b) Triangulation 
of the deformed atomistic configuration to visualize the Galerkin projection 
described in £2.3.1 The positions of the large red atoms are free, while the 
positions of the small black atoms are constraint by the motion of the free 
atoms. (c) Visualization of the QCE method described in £2.3.2 The blue 



shaded region is the set Vli , the red atoms inside the white region are the set 
^a Ce (both after deformation). (d) Visualization of an interface correction 



as described in §2.3.4 The blue shaded regions is the continuum region f2 c , 
the green shaded region is the interface region fii, and the white region is the 
atomistic region Sl a . The mesh in U fi a is chosen so that it coincides with ^ 
(cf. gO. 



For future reference, let 1^ : <3f — > &h denote the nodal interpolation operator. We note 
that 4 : W A ->• as well as I h : -> We also define Pf (5^) to be the set of all 
2Z d -periodic functions v h G Pi(^ # ); i.e., % = Pf {,%) d . Let &f denote the set of closed 
edges of the extended triangulation 5^ , and let denote the set of all edges / G J^/f such 
that area(/ n O) ^ 0. Finally, we introduce the spaces of piecewise constant functions Po(=^), 
Po(<5t ), and Pq (^), defined in a similar manner. 

The Galerkin approximation of ([7]) is the coarse-grained minimization problem 

y aj/l Gargmm< ot (y fe ), (10) 

Vh^Ath 

where we recall that <f a tot = S & + 
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For example, if we choose the triangulation 3~h as in Figure |ljb), then we obtain full atomistic 
resolution in the neighbourhood of the defect, while considerably reducing the overall number 
of degrees of freedom by coarsening the mesh away from the defect. In this way it is possible to 
obtain highly accurate approximations to nontrivial atomistic configurations. Under suitable 
technical assumptions it is not too difficult to employ the classical techniques of finite element 
error analysis in this context. Such analyses, including a 'posteriori grid generation, are given 
in 



Remark 1 (Regularity of Atomistic Solutions). In order for the Galerkin projection 
method, or the subsequent a/c coupling methods, to be accurate we require "regularity" of 
atomistic solutions. Such a regularity theory does not exist at present, however, most nu- 
merical experiments performed on atomistic models for simple lattices indicate "smoothness" 
of atomistic deformations away from defects. The situation would be different for so-called 
multi-lattices, which require a homogenisation step and represent a far greater challenge. □ 



2.3.2. Continuum region & Cauchy-Born approximation. The Galerkin projection (10) reduces 
the number of degrees of freedom considerably, however, the complexity of computing ^ a \w h is 
not reduced in the same manner. Due to the non-locality of the atomistic interaction ^ a |#h 
cannot be evaluated as easily as in the case of finite element methods for continuum mechanics. 



Several attempts have been made to use quadrature ideas to approximate S & and render ( 10 ) 
computationally efficient [18 } I21 | |2~I]. however, it was shown in |31j that these approximations 
yield unacceptable consistency errors. 

An alternative approach, proposed in [37j, is to keep the full atomistic description for atom- 
istically fine elements, while employing the Cauchy-Born approximation for coarser elements 
as well as an interface region. Following the terminology of [12] we call the resulting method 
the QCE method (the original energy-based quasicontinuum method). 

To formulate this method we choose a set 5£& cc C =5f of atoms that we wish to treat atom- 
istically (the red atoms in Figure [l|c)). Let Q £ (x) := x + e(— |, and define 

and n^ cc :=n\ [j Q £ (x). 

With this notation, the QCE energy functional is defined, for y^ G as 

GvxiVh) ■= e d V V{D m y h (x))+! W(Vy h )dx, (11) 

where W : ~R dxd — > M. is the Cauchy-Born stored energy function, 

W(F):=V(F<%) = V((Fr) re &). (12) 

Note that W^(F) is the energy of a single atom in the Bravais lattice FJzf*. 

One may readily check that the complexity of evaluating <# qcc (or its derivatives) is of the 
order 0(#^), that is, of the same order as the number of degrees of freedom. Moreover, it 
is easy to see that ^ qC e(?/A) = ^(ua) for all A G ~R dxd . For future reference, we give a formal 
definition of this property. 

Definition 1 (Global Energy Consistency). We say that an energy functional $ G C(^) 
is globally energy consistent if 

<% A ) = VA G R dxd . (13) 
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More generally one can show that ^q Ce (y/i) is a good approximation to # a (y/i) if ^Vh varies 
only moderately. Despite these facts, it turns out, as we discuss in { 2.3.3| and §4.2[ that 
minimizers of <# qce are poor approximations to minimizers of 

Remark 2 (The Cauchy Born Model). The Cauchy-Born model is a standard con- 
tinuum model, for large deformations of single crystals. In the absence of defects, solutions 
of a pure Cauchy-Born model (no atomistic region) can provide excellent approximations to 
solutions of the atomistic model ([7]). For example, in [17] it is shown that, for smooth and small 
dead load external forces, and under realistic stability assumptions on <f a , there exist solutions 
y a of ([7]) and y c of the Cauchy-Born model, such that 

( ed E E \ D *M X ) - D ej y c {x)A 1 ' 2 < Ce\ 
^ iey j=i ' 

where C depends on higher partial derivatives of V and on the regularity of y c . Hence, in 
the regime of "smooth elastic" deformations, the Cauchy-Born model can be considered an 
excellent approximation to the atomistic model ([3]). □ 

2.3.3. The patch test. The patch test is often employed in the theory of finite element methods 
[53\ O [5] as a simple test for consistency. The test also plays an important role in the design 
of a/c methods. 

Definition 2 (Patch Test Consistency). We say that an energy functional $ G C 1 (%;M) 
is patch test consistent if it satisfies 

(5<?(i/F),u fc )=0 Vn^tft, VF£l dxd (14) 
The terminology "patch test consistency" is motivated by Proposition [TJ where we have 



shown that the exact energy <f a does satisfy the patch test (14). 

However, the QCE energy (oqce is not patch test consistent |50j . This result will be reviewed 
for a one-dimensional model problem in §4.2[ where it will also be shown how failure of the 
patch test affects the consistency error. 

Remark 3. In most of the a/c coupling literature the patch test is stated as the condition 
that 

dyhip) 

It is straightforward to see that this condition is equivalent to the variational formulation given 



= for all finite element nodes p. 

Vh=yF 



in (14). □ 



2.3.4. Interface correction. In the engineering literature (see, e.g., |50| I35j) the non-zero forces 
under homogeneous strains of patch test inconsistent a/c energies are usually dubbed "ghost 
forces" . The discovery that <£q Ce is not patch test consistent has resulted in a number of works 
constructing new a/c methods that removed or reduced the "ghost forces" [4"1 \15 \ l4*9l [52"1 [56l |23"1 
[27 1 155 1 f22] . In some cases, this is achieved through sacrificing a variational (i.e., conservative, 
or, energy-based) formulation [TOl [201 GSl [331 S3 EH EE] . 

In the present paper we will focus only on energy-based a/c methods that are patch test 
consistent, i.e., that remove the "ghost forces" altogether. None of the methods presently 
available in the literature have resolved this problem in its full generality, however, |151 l4*9l [521 
[22] present several interesting approaches and partial solutions. In the following, we present a 
generalization of the geometrically consistent coupling method [15], which is the most general 
approach but leaves some questions concerning its practical construction open. 
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Let C 3?h and ft c = U^ c . We assume that all atoms belonging to Jzf \int(O c ) are vertices 
of £?h, and we define two sets J2? a c , ^ gc such that 

J^ gc UJ2f c = jSf \ int(n c ), and ,^ gc n J£f c = 0. 

For each x G Jz?j gc , we define a modified interaction potential V(x; •) G C 2 ((M d )'^), and we 



define the a/c energy functional as 

,(y h ):=e d £ y(Z^y(x))+£ d y(x;Lky(x)) + / W(Vyfc) dx. 



(15) 



For the modified potential V one takes a general ansatz with several free parameters, which 
are then fitted to remove or minimize the ghost force. For example, following the ideas of the 
quasinonlocal coupling method [52] and the geometrically consistent coupling method |15j one 
may define 

V(x;g) = V((g r ) r ese), where g r = ^ C x , r)S g s . 

The constants C Xt r )S can then be determined analytically as in |15[ I44j . or, as proposed in [45] , 
numerically in a preprocessing step. The 2D numerical experiments performed in (35] suggest 
that it is always possible to determine parameters C Xtr ^ s such that <^ ac becomes patch test 
consistent, however, a proof of this fact is still missing. 

The purpose of the present work is to investigate the question whether patch test consistency 
is in fact a sufficient condition for first-order consistency of an a/c coupling method. If this 
would turn out to be false in general, then it would be necessary to develop new approaches 
for constructing accurate a/c methods. 

2.3.5. General assumptions on the interface correction. For the subsequent analysis we assume 
an even more general form of the a/c functional than (15). We choose ^ c , S?£ C 
mutually disjoint, such that = U U J^ a , and we define the continuum, interface, and 
atomistic regions 

n c := U5^ c , 0; := U5£, and ft a := U5^ a . 

(Note that f2 c , f2j, J7 a are closed sets.) Next, we define the set of all nodes _£? a C J*? that interact 
with the atomistic region: 



J£f a := x G ££ 



z,x + er) n Qf / for some r G ^ j , 



where the ordered pair (x, x') G x Jzf# is called a bond; here, and throughout, the symbol 
(x,x') is also identified with the segment conv{x,x'} (in particular, it is closed). To avoid 
interaction between J2f a and Q c , we assume throughout that 

{x + tr \x G JS? a ,t G [0,1], r G ^} C ft a U fti. (16) 

Next, we define the set of interface bonds 

St- X := | b = (x,x + er) x G J??, r G ^, (x,x + er) C ft*}. 

Finally, we define an interface functional S\ G C 2 (^) such that 

^i(y) = e d E i ((D r y(x);(x,x + er) G #i)Y (17) 



that is, the interface functional <§\ is given as a function of the finite differences D r y(x) of bonds 
(x, x + er) that are contained in the interface region. Note also the volumetric scaling e d . 
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With this notation, we set 

^ac(lto) := e d V(D M y h (x)) + / W(Vy h ) dx + %(y h ). (18) 



The interface functional Si{yh) specifies the different variants of a/c methods. It is easy to 
see that the functionals ^ qce and <& c , discussed above, fit this framework (in the case of <^ qce 
we have to drop the assumption ([To])). 

If we define the total a/c energy as <f a ot := <^ac + ^aci where =^ ac is a suitable a/c approxi- 
mation to , the a/c approximation to ([7]) is 

y ac G argmm<4°% fo ). (19) 



If y ac solves (19), then it is a critical point of <# a ° = ^ ac + S ^: 

<^ ac (yac) + 5^ ac (y ac ), n h > = Vu h G <2T h . (20) 

2.3.6. T/ie locality and scaling conditions. We define notation for first and second partial deriva- 
tives of E\ as follows: for g = (gb)be&i let 

dbEi[(g b ) beS g.) := S and d b dy Ei^)^.) ■- 



lb dg b dg b > 

We extend the definition periodically: if b G ^ and £ G 2Z rf then d^ b Ei := d b Ei, and we make 
a similar definition for the second partial derivatives. In our analysis we will require two crucial 
properties on E\, which we call the locality and scaling conditions: 
The locality condition 

a r\ f \ n for all bonds (x,x + sr), (x' ,x' + es) G S$\ / 01 x 

d{ X , x +er)<){ X >, X '+es)E-Ay) - ^ ^ ^ ^ (21) 

implies that the same bonds interact through S\ as in the atomistic model. This condition can 
be weakened, by requiring that only bonds within an 0(e) distance interact, however, such a 
more general condition would add additional notational complexity. 

In the scaling condition we assume that there exist constants M* s > 0, r, s G St, such that 

( M* a , V(x,x + er),(x,x + es) £ 
|| WO < I ^ . f length( ^ f n (xj x + £r)) > , (22) 

This condition effectively yields an O(l) Lipschitz bound for 5S\ in the function spaces we 
will use. The scaling aspect enters through an implicit assumption on the magnitude of the 
constants M* s , namely, we will assume throughout that the constant 

Mi: =EEl r IN M ^ (23) 

is of the same order of magnitude as the constant M a defined in Q. 

Remark 4- The factor i for bonds on the boundary of the interface region is not strictly 
necessary, since it can be removed by simply replacing M* s with 2M l r s , however, if stated as 
above it makes the statements of the results in ^slightly sharper, and moreover simplifies the 
argument in (68). 

The necessity of this factor is related to the fact that we allow E- x to depend on bonds that 
lie on the boundary of fij ; this is made clear in Proposition 
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where we construct a stress 
function for $ &c . Note that if we did not allow E\ to depend on these boundary bonds, then 
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it would in fact be impossible to construct patch test consistent a/c methods for non-flat a/c 
interfaces. □ 

3. A Framework for the a priori Error Analysis of a/c Methods 

When analyzing the error of a numerical method, one should first of all determine the main 
quantities of interest. For a/c methods, one is usually interested in energy differences between 
homogeneous lattices and lattices with defects, or critical loads at which defects form or move 
(i.e., bifurcation points). Since the present paper is mostly theoretical, we will simply focus on 
the error in the deformation gradient. We note, however, that many aspects of this analysis 
are crucial ingredients for the analysis of energy differences (see, e.g., [H] ) and would usually 
also enter an analysis of bifurcation points. 

We assume from now on that d S {1,2}. To execute the abstract framework of this section 
also in 3D, several technical tools as well as the central consistency result need to be developed 
first. 

3.1. Discrete and continuous functions. In the following analysis it will be important to 
extend the a/c functional <# ac to all functions y E & '. To that end, we first define piecewise 
affine interpolants with respect to an atomistic mesh 3? £ . 

We take a subdivision of the scaled unit cube e(0, l) d into d-simplices (in ID the interval 
e(0, 1); in 2D two symmetric triangles; compare with the triangulation of the atomistic region 
in Figure [l|d)), which we extend periodically to a triangulation 2?^ of M. 2 with vertex set 
«Sf*. The restriction of to $7 is denoted by 8F e . Each discrete function v : — > R fe will 
from now on be identified in a canonical way with its continuous piecewise affine interpolant 

v e Pi(^ e # ) fe . 

For future reference w e deno te the sets of edges of 2? e and 2?^ , corresponding to the defini- 



tions of & h and in 5 2.3.1 by and &l 



3.1.1. Ambiguity of continuous interpolants. If € then can also be interpreted as a 
member of *3f and therefore has two, possibly different, continuous interpolants. To distinguish 
them, we make the convention that the symbol yh always denotes the interpolant in Pi(^) d , 
while a symbol y always denotes the interpolant in P\(£7 £ ) d . If we wish to evaluate the Pi(£7 £ ) d - 
interpolant of a function jfe e % then we will write I e yh- 

To compare a Pi (^^-interpolant with a Pi(^) d -interpolant, we use the following lemma. 
In ID the result is easy to establish; in 2D it depends on a technical tool that we introduce in 



£5.1 The proof is given in the appendix. 



Lemma 2. Let d £ {1, 2}; then, for all yh £ &h and p G [1, oo], we have 

||VI £ yft||Lp(n) < Wyh\\w(Q)i 
where we recall that we have defined ||Vu||lp = || |Vu| p ||lp- 

3.1.2. Extension of the a/c energy. Before we extend the a/c energy <f ac we make one last 
technical assumption, which considerably simplifies the subsequent analysis. We shall assume 
from now on, that 

If all atoms in ££ n (O a U Oj) are vertices of ^ a U which is not uncommon, then this is no 
restriction. 
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With these conventions the a/c energy <# ac defined in (18) can be defined canonically for 
functions y = y + Uh £ & + % by the same formula: 

^ac(y) = I W(Vy) dx + e d J2 V{D M y(x)) + S^y) forjefi W h . 
We also assume that & ac has a suitable extension to It should be stressed that for general 

3.2. Measuring smoothness; An interpolation error estimate. The three main ingre- 
dients in the a priori error analysis of Galerkin-like approximations are (i) consistency, (ii) 
stability, and (hi) an interpolation error estimate. We begin by establishing the latter. To that 
end, we first need to find a convenient measure of smoothness for discrete functions y E W . 

3.2.1. Measuring smoothness in terms of local oscillation. There are several possibilities to 
measure the "smoothness" of a discrete function. The most obvious is possibly the use of higher 



order finite differences, e.g., D ei D e .y(x). If, in £3.1 we had chosen continuous interpolants 



belonging to W ,0 °, then we would be able to simply evaluate the second derivatives V y. 
However, since the interpolants we use are piecewise affine, the second derivative of y is the 
measure [Vy] ® i/ds\^.#, where [Vy] denotes the jump of Vy across an element edge, and ds 
the surface measure. 

This last observation motivates the idea to measure smoothness of y by the local oscillation 
of Vy. We define the oscillation operator, for measurable sets ui C M. d , and for y S & , as 

(x? ^ \Vy{x)-Vy(x')\ 
osc(Vy; u) := ess sup -. (24) 

x,x'£u) ^ 

The sets u> that arise naturally in our analysis will always have 0{e) diameter, which is the 
reason for the e _1 -scaling in the definition of osc. 

Note, in particular, that if y were twice differentiable, and if diam(w) < Ce, then we would 
obtain 

diam(u;) 2 

osc(Vy;a;) < ||V y\\L°°(u) < C\\V y|| L °°(w), 

which further illustrates that the oscillation operator is a reasonable replacement for V 2 y to 
measure the local smoothness of a piecewise affine function. 



3.2.2. Interpolation error estimate. The smoothness measure we defined in §3.2.1 yields a sim- 
ple proof of an interpolation error estimate; see Appendix [A) 

Lemma 3. Let d G {1,2} and suppose that ^ U Sft C S? e ; then there exists a constant Ci, 
which depends only on the shape regularity of such that, for all y E W , p S [l,oo), 



l/p 



f r ip} 

||V(y-4y)|| LP(n) <CA Yl |r|[h T osc(Vy;^)J \ 

where, for T 6 £? £ , 

^■.= nfn\J{T' e ST*\TV\T' ^0}, and h T :=max|/i(x)|. (25) 
Similarly, for p = oo, we have ||V(y — ^i2/)|| L <x>(n) < C/max-pg^c [h-p osc(Vy; u>^)] . 
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3.3. The stability assumption. The stability of a/c methods relies, firstly, on the stability 
of atomistic models as well as their Cauchy-Born approximations. It requires a thorough 
understanding of the physics of a model and in particular more specific information about the 
interaction potential. Since the focus of the present work is the consistency of a/c methods, we 
will formulate stability as an assumption. 

The simplest notion of stability one may use, which is also closely connected to local mini- 
mality, is coercivity of the second variation: 



(5 2 4°\y h )u h ,u h ) > co\\Vu h 



2 



V«h E & h , 



(26) 



where cq > 0, and yh E % is a suitable deformation in a neighbourhood of the atomistic 
solution y a , e.g., yh = IhU&- The choice of norm is motivated by the fact that the Cauchy-Born 
model, and hence the atomistic model, are closely related to second order elliptic differential 
equations. 

Examples of sharp stability estimates for a/c methods in ID can be found in |13 } 1401 1^71 I26| 
I55j . For pair interactions in 2D the stability of Shapeev's method [49J is established in |41j . 

More generally, for some p E [l,oo],p' = p/(p— 1), we may assume an inf-sup condition of 
the form 



inf 

u h e<& h 

|[V«ft||LP(tJ) = 



sup 

:1 l|v^|| L y (n) = 



(6 2 <?£\y h )u h ,v h ) >c , 



(27) 



for some constant cq > 0. The condition (27) is usually difficult to prove, especially for p ^ 2, 



and may indeed be false in general. We will only use it to demonstrate how such a stability 
result motivates consistency estimates in different negative norms. Examples of ID inf-sup 
stability estimates for a/c methods can be found in [TO l [T4"l I 



3.4. Outline of an a priori error analysis. The following outline of an a priori error analysis 
depends on a stability assumption that we will not prove. Moreover, since it primarily serves 
to motivate the consistency problem, and since a rigorous derivation would be more involved 
without yielding much additional insight, some steps will be kept vague. Most of these steps 
are easily made rigorous; the main assumption we make below, which is in fact very difficult 
to justify rigorously, is that lhV& and y ac are "sufficiently close". See [32l HQl H21 [55] for similar 
analyses in ID where all steps are rigorously justified, and [43^ |4"T] for a similar semi-rigorous 
framework, where a proof of this step is replaced by an assumption. 
Let y a satisfy ([8 
I h y a . Let e h : 



y ac satisfy (|20|), and suppose that the stability assumption (|27|) holds with 



Vh = ihVs.- -Let e h := y ac - i h y&- Moreover, suppose that \\VI h y, 
so that the following approximation can be made precise: 



Vy a 



L°° 



is sufficiently small 



5 2 S^{I hyi , + te h )e h ,v h )dt 



(5(<f ac + ^ ac )(2/ ac ) - <K^ac + ^)(hy a ),v h ). 



Taking the supremum over all Vh E and invoking the inf-sup condition (27) and the criti- 



cality condition (20), we obtain 



where we define 



< |K^ac + ^ac)(4ya)|| w -i,P, 

|$|| w -i, P := sup ($,«/,), for$€%*. 



\\Vv h \\ LP ,=l 
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We split the consistency error ||<5(<£a C + ^ 2l ac)(Aya)||-w~ 1 ' p t nree separate contributions: 



h 



^11^11^(0) ~ W 5 ^ c + ^*ac)(/h!/a) " + ^ac)(ya)|| w -i,P 

+ ||^ac(ya) " ^a(ya)|| w -l,P + ||<^ac(ya) " <S£»a(l/a) || w -i, P 
. c*coarsc , cmodel , cey± 

where we have used (Tsl) , and the extension of <^ ac and for all deformations y G ^ constructed 



in S 3.1 



The coarsening error, £^ oarse ; can be bounded by Lipschitz estimates for 5(<£ac + <^ac) and 
an interpolation error estimate. Using our assumption that U C =5^ it is not difficult to 
derive Lipschitz estimates of the form 

£ coarse < ^ + M ^j|| V I ft y a _ Vj/ a 1 1 ^ , (29) 



where M a is a Lipschitz constant for dW (cf. (12)), and M^ ac is a Lipschitz constant for <5^ ac . 
Combining (28), and (29), the inequality 

||Vy a - v yac|| LP < ||Vy a - V4y a || LP + ||Ve fe || LP , 

and the interpolation error estimate of Lemma [3| we arrive at the following basic error estimate 



i/p 



(30) 



l|Vya-Vy ac || LP(n) < h J h +-U £ [h T osc(V ya ;^)] 

where a = Cj(c + M a + M^ ac ). 

The consistency error for the external forces, £^ xt , depends on the form of & a and <$^ ac and 
cannot be discussed at this level of abstraction. The modelling error, 

£ mode\^ is the focus of the 

remainder of the present paper. 

Remark 5 (Choice of Splitting). Suppose, for simpliciy, that <5^ ac = & & = 0. In 
a typical finite element error analysis of continuum mechanics problems one would usually 
choose a different splitting of the consistency error: 

||^ac(4ya)|| w -i,P < ||^ac(4y a ) ~ ^a(42/a) || w ~hP + ||^a(4y a ) ~ 5<f a (y a ) || w -i,p. 

This splitting was used in the analysis in [41] and led to a suboptimal estimate of the modelling 
error, since it still contains some coarsening error. □ 

3.5. The consistency problem. The main step that remains in obtaining an a priori error 
estimate from (30) is the estimation of the modelling error 



£h = *<Miw - ^a(y a ) w -i,p = sup -— . 

h v h eW h \{0} ll v ^llLp'(n) 



Most of the numerical analysis literature on a/c methods estimates this modelling error only 
for the case when ^ = 3T e . In ID it is easy to see that this is sufficient, since I £ Vh = Vh in that 
case; see also [33]. The following lemma provides the main technical step to explain why it is 
also sufficient in 2D to consider the case ^ = Its proof uses arguments similar to those in 
the a posteriori error analysis of continuum finite element methods and is given in Appendix 

Lemma 4. Assume that 5^ U S?£ C Let $ G ^* and $^ G ^* be given in the form 

= / a : Vudx, and ($^,1%) = I a:Vuhdx, 
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for all u G % u/t G where a G P*(^) rfx<i ; i/ien £/iere exists a universal constant Cm such 
that, for all p G [1, oo), 

/ \ i/p 

($,I e u h ) -{§h, Uh) <C M e( ^ \T\osc(a; ujt) p j \\^u h \\ LP ,^y and 

^Te^ ' (31) 



%I e u h ) - (<S> h ,u h ) 



< 



Cm? \ max oscfcr ; VitJ T1 



(He)" 



The estimate (31), together with Lemma [2j implies the following theorem, where we use the 
notation 



l$l 



sup (3>, 

ve<&\{0} 
l|V«|| LP ,=l 



for 3> G <2r*. 



(32) 



Theorem 5. Suppose that 2?£ U J,' C 5 £ and f/iai y 6 ®0 ^en, /or aZZp G [l,oo) ; 

||<5<? ac (y) - 6^(y)\\ w - hP < M a C M s( ^ |T|osc(Vy(T); oj^)A * 



(33) 



Te^ c 

+ ||^ac(y) - <^a(y)|| wr i, P , 

lyii/t corresponding statement for p = oo. 

Proof. Since S a {I e yh) uses only point values of I e y^,, which are the same as for y^, we have 

(<5^ a (y),«fc> = {S^(y),I £ u h ) V« A G 
Using this fact, we can estimate 

|(^ac(y) - fi&a.(y),Uh)\ < |<^ac(y),^} - (8&ac(y),IeUh)\ + \($&ac(y) ~ 5£ a (y), I £ U h )\. 

Due to the assumption that 2?£ U C the first group can be estimated using Lemma |4| 
with a = dW(Vy), which yields the first term in (33). 
Using Lemma [2j the second group can be estimated by 

\(S&ac(y) - S^a,(y),I e Uh)\ < \\S&ac(y) - S&a.(y)\\ W jl,p\\VI e Uh\\ IJ / 

< ||^ac(y) - s ^(y)\\ w -i-p\\^ u h\\ LP > ■ 

Taking the supremum over all Uh G % with ||V«/j,||t p ' = 1 yields the stated result. 
Applying Theorem [5] to the modelling error £™ odel , defined in (28), we obtain that 



□ 



^modci < f modd + M a C M e{ \T\osc(Vy(T);uj^) p 



i/p 



(34) 



where 

C° dCl :=||^ac(ya)-^a(ya)|| w -l,P. 

Even though £™ odel - g e g sen tially an upper bound for £™ odel , it is usually easier to estimate. 
The consistency problem is to prove a sharp upper bound on £ ™ odel . 

In £j4] we will discuss two simple ID examples to determine what can be expected in more 
general situations. In Theorem 10 we will prove that for an a/c method that is patch test 
consistent, and satisfies various other technical conditions, one obtains 

i/p 



£model ^ 



f 1 /P 

Cel Yl m[osc(Vy a ;^ T )n , 
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where C is a constant that is independent of y a , but does depend on the interface width, and 



ojt C r2 c U fii is the interaction neighbourhood defined in (43). 



Combined with (34) and (30), and using the fact that e < hy, and that uy D 10 r ^ £ <^£ C > 
this bound yields 



Vy a - Vy a 



LP 



< 



cxt 



— +^ 

CO c 



E 



ITI 



h T osc(Vy a ; u T ) 



l/p 



(35) 



where C2 is a constant that is independent of y a . This estimate closely resembles a typical 
first-order a priori error estimate for a continuum mechanics finite element approximation; see 
also the interpretation given in j |1.2| 

It should be stressed again that (35) is not a rigorous error estimate, but depends on various 



assumptions made in the forgoing discussion, most prominently, the stability assumption (27), 
and the assumption that \\VIhVa, — Vy ac ||L°° is "sufficiently small". 



Remark 6. The locality of the patches lot is crucial. If diam(wj') is not of the order 0(e), 
then it is possible that osc(Vy a ; wt) 3> 1 even if y a is globally smooth; see also £6.4.4 □ 



4. Examples in ID 

In the present section we review the consistency analyses of specific a/c methods to point 
out the main features and to motivate what may be expected in the general case. Throughout 
this section we assume that d = 1, ffl = {±1, ±2}, and that V is given by 

V{{g±i,g± 2 }) = §[0iOi) + Md-i) + M92) + ^2(9-2)], 

where 4>i,4>2 £ C 2,1 (1R) are, respectively, the first and second neighour interaction potentials, 
which are assumed to be symmetric about the origin. We assume that their derivatives (\)\ and 
(j)'l have global Lipschitz constants m! i and m". 

For the ID analysis it is convenient to write x n = ne, v n = v(x n ), and to write all interactions 
in terms of the backward difference operator 

/ v n - V n -l 



£ 



With this notation the atomistic energy can now be rewritten in the form 

Af N 

<^a(y) = e My' n ) + z Yl My'n + y'n+i), 



(36) 



n=-N+l 



n=-N+l 



where we note that y' n + y' n+ i = £ 1 (|/n+i — 2M-i) describes a second neighbour bond. 
For future reference we also define the second and third finite differences 



J n+1 



and v' n '(x) 



V 'n+1 ~ 2v 'n + v 'n-l 



It is also worth pointing out that v' n = Vv(s) for all s G (x n -\, x n ). 

4.1. Consistency of the QNL method. We begin with a modelling error analysis for the 
quasinonlocal coupling method (QNL method) of Shimokawa et al [52j. The geometrically 
consistent coupling scheme [15] and the method proposed by Shapeev [19] reduce to the same 
method for ID second neighbour interactions. 

The following presentation follows largely [40] . where the QNL method is defined as follows: 
Let Ji = {-K, . . . , K} for some K > 1 and jV c = {—N + 1, . . . , N} then, for y6f, the 
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QNL energy is defined by 

N 

^qnl(y) = e MVn) + £ £ ^« + ^n+l) + £ E 2 MV„) + M^'n+l)} ■ (37) 

We observe that we have not modified the first neighbour interactions, but have "split" the 
non-local second neighbour interactions into local first neighbour interactions in the continuum 
region. 



It is straightforward to rewrite ^ qn i in the form specified in (18), with 



n c = [e(K + l),l}U(-l,e(-K-l)}, jg? a = e{-K + 1, . . . , K - 1}, 



and a suitably defined interface functional $\; however, the form (37) is more convenient for 
the analysis. 

The following modelling error estimate was first established in |40[ Thm. 3.1]. Dobson and 
Luskin [12] treated a quadratic interaction case, using entirely different analytical techniques 
that gave an even more detailed analysis of the error; Ming and Yang [36] used related methods 
as |40| Thm. 3.1], but gave a qualitatively less precise estimate of consistency error. An 
extension of the result to linear finite range pair interactions is given in [27] . 

Note also that it is shown in [12] that the consistency error of the QNL method in f e -norms 
is of the order 0(1), that is, the usage of negative norms cannot be avoided. 

We will discuss the estimate in detail in §4.3[ The proof of the following result serves as a 
first guidance on how one may approach proofs of consistency of a/c methods in more general 
situations. 

Proposition 6 (Consistency of the QNL Method). Let y 6 f ; then 

||^qni(y) - &&(y)|| w -i,P < £ml 2\\y"\\p e{{ ^ K)K }) + ^^WW^^i) + ^^Wy'^p^y 

where = {—N + 1, . . . , —K — 1} U {K + 2, . . . , N}, 

Proof. Throughout the proof we will make use of the fact that the boundary conditions are 
periodic without comment, treating the boundary as if it belonged to the "interior" of the 
continuum region. 

Since the first neighbour interactions as well as the second neighbour interactions in the 
atomistic region are treated identically in the atomistic model and the QNL method, we have 

(Vn + Vn+l) ' ( u 'n + u 'n+l) 



^(Zy'n) ■ u 'n ~ ^y'n+l) ' «n+l 



Rearranging the sum in terms of the gradients u' n , and using the fact that Jf^ = {— K, . . . , K}, 
yields 

N 

(S^(y)-5^ qnl (y),u h ) = e £ R n ■ u' n , (38) 

n=-N+l 

where (R n )^ =1 is defined as follows: 

0, ne{-K+l,...,K}, 

Wn-l + Vn) ~ M 2 l/n)> n = K + 1, 
Mn+y'n+J-tiPy'n), n = -K, 
{ Wn + y'n+l) ~ WWn) + Wn-1 + t/nl n * ■ 

At the interface, n = K + 1, we have 

Rn < m' 2 \y' n+1 - y' n \ = em'J^J, 



R> 
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with a similar estimate for n = —K. In the continuum region, n > K + 2, or n < —K + 1, the 
terms R n have second order structure, and a second order Taylor expansion yields 



|Rn| < \4>2{2y l n )\\y' l -2y' n + y' n _ 1 \ + 



l // 

2 m 2 



\ i / I i I i f 

\Un+l ~ Vn\ ~+~ \Vn ~ Vn x \ 



2 1™// 



£ m 2\Vn \ + e 2 m 2 



Wr. 



+ y: 



n-l 



After inserting the two bounds into ( 38 ) , and applying several weighted Holder inequalities one 
obtains the stated estimate. □ 

4.2. Inconsistency of the QCE method. As in the previous section, let ,JV & = {—A, . . . , A}, 
Jf c = {— N + 1, . . . , N} \ and J2? a = e,jV & \ then, under the assumptions and notation set out 
at the beginning of the QCE energy functional defined in (11) reads 

^qce(y) = E ^2 5 [MVn) + MVn+l) + fail/n-l + Vn) + fail/n+1 + ^+2) 



/ / 2 



(39) 



noting that W(F) = <M F ) + <fo( 2F )- 

The following result is a variant of [431 Thm. 3.2]. Previous analyses of the QCE method 
Q21 [36] computed the consistency error contribution due to the "ghost forces" explicitly 
rather than estimating their ^-residual contribution. 

Proposition 7. Let y G <3f and p G [l,oo]; f/ien 

||^qce(y) - 8&*(y)\\ w -i,P < el/PG + £m 2\\y"\\ £ P^) 



i 2 / II /// 

+ e m 2 2/ 



I 2 //II //|| 2 P 

m 2||y ||^(xu{-jc,jc})' 



where j¥£ is defined as in Proposition^ J\f{ = {—A — 1, —A, K,K + 1}, and 
G = |[|^ 2 (2y^_ 1 )| + \4>W-k + i)\ + \<&W K )\ + \<t>' 2 (2y' K+2 )\]. 



The estimate is sharp in the sense that, for some constant C , \ < C < 2, 



|^qce(y. 



qceli/Aj ||\y _1 ' p 



W, 



> 



Ce 1 /P|^ 2 (2A)| 



VA G 



Proof. The first result can be proven in much the same way as Proposition |6j by rewriting the 
first variation 5S qce in the form (&f qce (?/), u) = e Y^, n =-N+i ^™ ' u 'm an d carefully estimating the 
coefficients R n . See |43j for the details of this computation. 

To obtain the opposite estimate for y = y/\, a brief computation gives 

(5^qce(?/A) - &f a (2/A), = £2^2(2A) ~ U_ K+1 ~ u' K + u' K+2 ] . 

If we choose u£f such that 

( (4e)- 1 / p ', n = -K -l,K + 2, 

< = sign(0 2 (2A)) • ^ -(4e)- 1 /p', n = —A + 1, A, 

( 0, otherwise, 

then ||Vu|| LP ^ _ 1 ^ = ||n'||^ p ' = 1, and we obtain that 

H^qcc(yA) " <S&(!/A)|| W -l,i. > <^qce(y A ) " 6£ a (y A ), U ) = 6*'* [2 ■ ^ 1 0' 2 (2 A) |] . □ 
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4.3. Discussion. This discussion of the ID consistency error estimates largely follows the 
discussions in [131 00] . 

We have estimated the modelling errors for two prototypical a/c methods. We see that the 
leading order terms in the upper bounds are O(e) and 0{e l / p ) for the QNL and QCE methods, 
respectively. However, a much finer distinction should be made. 

First, we note that both methods reduce to the Cauchy-Born approximation in the continuum 
region, and the corresponding contributions are all of second order (see also Remark [2j note 
also that this requires point symmetry of V, which we have not assumed in general in this 
paper). 

Second, we see that the QCE method (and only the QCE method) has a zeroth-order term 
Ge l l p in the interface region. This term occurs since the QCE method is not patch test 
consistent, that is, homogeneous deformations are not equilibria of the QCE model: 

^qce(yA) / 0. 

The origin and effect of these "ghost forces" are discussed in more detail in [5Ql ESI Ell 021 E6j • 
We should call this term zeroth order for several reasons: Firstly, it is clearly of zeroth order 
lip = oo, in which case the consistency error is related to the error in the W 1,00 -norm. Secondly, 
the parameter e is a constant of the problem and does not tend to zero. As a matter of fact, 
the accuracy of an a/c method should be related the smoothness of the solution (as opposed to 
the atomistic scale), and the term Ge 1 ^ is independent of the magnitude of y" in the interface 
region. The scaling e 1/,p relates only to the width of the interface region. 

Finally, it is worth remarking on the first-order consistency term in the interface region for 
the QNL method. The reason this term is of first order as opposed to second order is the loss 
of symmetry that is introduced by changing the interaction law at the interface between the 
atomistic and continuum regions. A recent result of Dobson [H] shows that no a/c method 
coupling an interatomic potential to the Cauchy-Born continuum model can achieve better 
than first-order accuracy in the interface region. 

Note also, that the second finite differences can in fact be written in terms of the oscillation 
operator: 

\y"\ = osc(Vy; [x n -e,x n + e]). 

In our analysis in we will ignore the possibility of proving a modelling error estimate that 
is of second order in the continuum region, but we will be satisfied with an estimate that is 
globally of first order. 



5. Auxiliary Results 

5.1. The bond density lemma. The bond density lemma is a tool that allows a transition 
between integrals over bonds, and volume integrals. It was first derived in [49J for the construc- 
tion of patch test consistent a/c methods for pair potential interactions. Related asymptotic 
results were used previously in T-convergence analyses of atomistic models [3] ; the achievement 
of Shapeev [39] was to obtain a formula that is exact for any triangle. 

Before we formulate the result we introduce some notation. Let T C K 2 be a triangle with 
vertices belonging to =§f*. We define the characteristic function xt : M 2 — ?• IR by 



Xt{x) := lim 



\TnB t {x)\ 



t\0 \B t (x)\ 

where Bt(x) denotes the closed ball with centre x and radius t. Note that xt = 1 in int(T), 
and xt = 1/2 on the edges of T. The value of xt on the corners is not of importance. 
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Let x,x' E M 2 and let if be a function that is measurable on the line segment (x,x'); then 
we define the line integral, or bond integral, 



I'x' pi 

f <pdb:= cp((l -t)x + tx')dt. 

Jx Jt=0 



»2 



Lemma 8 (Bond Density Lemma ([M], Lemma 2)). Let T C K be a non- degenerate 
triangle with vertices belonging to Jzf* = el? , and let r 6 Z 2 ; then 

fx+er 



e 2 



Y, f XTdh = \T\. (40) 

As a first application of the bond density lemma we present a proof of Lemma [2] in the 
appendix. 

Remark 7. In the above form, the bond density lemma is false in 3D, which is one of the 
reasons why the present work is restricted to 2D. Moreover, the condition that the vertices of 
T belong to lattice sites is also necessary. This is related to the assumption that the vertices 
of belong to Jz?*, however, this is not crucial and could be removed with some additional 
work. □ 

5.2. Discrete divergence-free tensor fields in 2D. Our second auxiliary result concerns 
representations of discrete divergence-free Po-tensor fields. Following the construction given by 
Polthier and Preufi |48j . we will give a proof for the periodic setting. This proof also serves to 



motivate a crucial argument in \ 6.3.3 Since we will only use the atomistic finite element mesh 
2T e from now on, we will formulate everything in terms of this mesh. However, all results hold 
for general periodic triangulations. 

5.2.1. The Crouzeix-Raviart finite element space. The representation of discrete divergence- 
free tensor fi elds requir es the use of the non-conforming Crouzeix-Raviart finite element space. 



Recall from £3.1 



2.3.1 



the definition of the sets of edges J£" £ and ■ The Crouzeix-Raviart 
finite element space over is defined as 

Ni(^*) = {w : U Tg ^#int(T) — > R | w is piecewise affine w.r.t. and 

continuous in edge midpoints ^t}- 

jj 

The degrees of freedom for functions w G Ni(^ ) are the values at edge midpoints, w{qf), 
f G i and the corresponding nodal basis functions are denoted by (f. 

The Crouzeix-Raviart finite element space of periodic functions is defined as 

Nf(5£) = {w e Ni(^ e # ) | w{£, + x) = w(x) for dx-a.e. x£R 2 ,(£ 21?}. 

jj 

The periodic nodal basis functions are defined, for / E J? e , by Cj = X^e2Z 2 

5.2.2. Path integrals. For two edges /, /' G ^fi let f denote the set of all piecewise affine 
paths from qj to qy, crossing element edges only in edge midpoints; see Figure [2| a) for an 
example. 

For any piecewise constant vector field a G Po(=^*) 2 and for any path 7 G ^fj', 7 = 
{x(t) I < t < T}, we denote the standard path integral by 

a ■ dx = a ■ x(t) dt. 

J-y Jt=0 
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Figure 2. (a) Illustration of a piecewise affine path 7 G Pf,f- (b) Illustration 
of the path j q used in the proof of Lemma [9j 

For piecewise constant tensor fields a G Po(=^*) fcx2 we define the path integral as 

a ■ dx = / ax(t) dt. 
7 Jt=o 

Since functions w G Ni(^*) fc have piecewise constant gradients Vw, and since they are 
continuous in edge midpoints, it is easy to see that 



/ Vw ■ dx = w(qfi) — w(qf) V7 G Ffj'- 
J 1 



(41) 



5.2.3. Discrete divergence-free tensor fields. The following lemma characterizes discrete divergence- 
free tensor fields. 

Lemma 9. A tensor field a G Pf (<%) kx2 satisfies 

[a:Vudx = Vu G Pf (,%) k 
Jn 

if and only if there exist a constant ctq G M fcx2 and a function w G N*(^) fc such that 



a = o"o + VwJ, where J 



-1 

1 



G SO(2). 



Sketch of the proof. The reverse direction, that any tensor field of the form a = o~o + Vu>J has 



zero discrete divergence, can be checked using a straightforward calculation, using (41). 

Step 0. Outline: To simplify the notation, we define a = aS T and, without loss of generality, 
we assume that k = 1. Initially, we treat a as a piecewise constant tensor field on all of R 2 , 
ignoring periodicity. We construct w G Ni(^*) by explicitly specifying w(qf) for all edges 

f G . We will then show in the last step of the proof that w can be written as the sum of 
an affine function and a periodic function. 

Step 1. Construction of w: Fix a starting edge / G ^ e and define w{q^) = 0. For any edge 

/ G let 7 G and define w(qj) via the path integral 



a • dx. 



We need to show that this definition is independent of the path. 
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Let q be a vertex of the triangulation ^ and let ip q be the corresponding nodal basis function 
with support uj q ; then a fairly straightforward calculation (see, e.g., [38] for the details) shows 
that 

0=/ a-Vcp q dx = - ^2 (( (T - u f) + + ( cr - ,/ f) ) = l a- Ax, 

q /e^# 79 

/Cint^q) 

where z^ 1 are the two unit normals to /, and 7 9 is the piecewise affine path through edge 
midpoints circling q; cf. Figure [2^b) . Note that the rotation J in the definition of a comes from 
the fact that tangent vectors are rotated normal vectors. 

Since all closed piecewise affine paths can be written as a sum over paths 7 g , this implies 
that J a ■ dx = for all closed piecewise affine paths 7, and in particular that the definition 
of w(qj) is independent of the choice of path, that is, w is well-defined. 

Step 2. V«) = a: From the definition of w(qf), f G , it follows that, for /, / C T G E? e , 

Vui(T) • (q f - q fl ) = w(q f ) - w(q f >) = a(T) ■ (q f - q f >), 

which immediately implies that Vw(T) = a(T). 

Step 3. Periodicity: We are only left to show that w(x) = a ■ x + w\(x) for some a G M 2 , 

W i{lf + 2ej) = wi(q^) = 0, we have 

w i(Qf + 2ej) = / a • dx = a ■ dx = wi(qj). 

J-y+2ej J 7 

This shows that w\ is periodic and thus concludes the proof. □ 

6. A General Consistency Result in 2D 

We are now finally in a position to make precise the statement that patch test consistent a/c 
methods are first-order consistent. Motivated by the example of the QNL consistency result (ig- 



a = (ai, 02), and wi G N^(^). Let aj = w(q^ + 2ej), j = 1, 2, and define wi(x) = w(x) — a ■ x. 
Fix j G {1, 2}, let / G and let 7 be a path from qt to qf. Since a is 2Z 2 -periodic, and since 



noring, as discussed in £4.3 the second-order consistency of the Cauchy-Born approximation), 
we would like to prove a result of the form 

\\5^{y) - <54(y)|| w -x, < Ce\\vM\ LP{ncUni) . 



As explained in £3.2.1 we will use the oscillation operator (24) to replace the undefined second 
derivative. 

For each element T G <5£ , let ojj, be the interaction neighbourhood of T in the atomistic 
model, 

uj* ■= [x + tin + t 2 r 2 I z G T, ti G [0, 1], r< G (42) 
and wt its union with the set w^, restricted to the continuum and interface regions, 



(43) 



TDTY0 



Note that assumption (16) implies that U C for all T G =^ c . 



Theorem 10 (First-order Consistency of Patch Test Consistent a/c Methods). 

Suppose that <# ac is patch test consistent ( §2.3.3 ) and globally energy consistent ( §2.3.2 ), that 
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the locality condition (21) and the scaling condition (22) hold, that int(O a ) is connected, that 
U J^ 1 C ZF e , and that (16) holds. Then, for any y G we /iawe 



<^ac(y)-^a(y)|| w -i, P <e\ Yl \T\\M T osc(Vy;uJT 



(44) 



where the oscillation measure osc is defined in (24), £/je interaction neighbourhood ux is defined 



in (43), and the pref actors Mt are defined as follows: 
M T = < 



0, T G =5^ a , 

(M i + M a )(l + 7width(O i )), Te^, 
M a + 7(M ; + M a )width(Oi), T G 



(45) 



In (45), the constants M a and M 1 are defined in §6§ and (23), and the "interface width" 
width(f2i) is given by 

length(^) 

width(f2j) := max min min . (46) 

/<=«# f cn f 

Outline of the proof. We will construct "stress functions" S a (y), £ ac (y) G Po(=^) 2x2 such that 

(5<f a (y),u} = / S a (y) : Vudx, and (<5<f ac (y); u) = / S ac (y) : Vudx Vu G . 
Jq Jn 

If we could prove an estimate of the form 

|£ a (y;T) - X ac (y;T)| < eosc(Vy;w T ), 

then the result would follow immediately. It turns out that this is not possible. 

Instead, we will use the fact that <f ac is globally energy consistent and patch test consistent 
to construct a correction (cf. Corollary 15, Lemma [TTJ and ^6.2.3 ) 



£ ac (y) = X ac (y) - W>(y)J, 

which still represents the first variation b<§ &c (cf. Lemma [9]), and which has the property that 

S a (y F ;T) = £ ac (y F ;T) = dW(F) VTG^U^, F G M 2x2 . 

In addition, we show that X ac (y; T) = £ a (y; T) for all T G 

Lipschitz estimates for S a and X ac , and careful modifications of the argument in the interface 



region, yield the following result (cf. Lemma 19): 

|S ac (y; T) - S a (y; T)| < |S ac (y; T) - dW{Vy{T))\ + |£ a (y; T) - dW(Vy(T))\ 
< eM T osc(Vy; uj t ) VT G & e , 

and, in particular, 

(S^(y) - ^a(y); u) = \ T \ \^c(y; T) - £ a (y; T)j : V«(T) 

i/p 



< el \T\ M T osc(Vy;ui T ) 

which yields the stated first order consistency estimate. 



Vul 



Lp' (fi) ) 



□ 



26 



C. ORTNER 



Remark 8. In ID, a similar result can be proven using a similar framework but with 
significantly reduced technicalities. Note, in particular, that the ID analogue of Lemma [9] is 

a ■ Vn dx = VuG^ if and only if a is constant. 



Hence, the corrector function ip(F, •) defined in £6.2.2 is always identically equal to zero, which 
removes the interface width dependence from the modelling error (cf. j ffi.3.3 ). Hence, if d = 1, 



we obtain (44) again but with modified prefactors 

f 0, TG^ a , 

A4 D := <^ M 1 + M a , T G (47) 
{ M a , T G 3? £ c . 

Moreover, since ip = in ID, and since symmetries are more easily exploited, it is not too 
difficult in ID to prove second order consistency in the continuum region. □ 

6.1. The atomistic stress function. A natural "weak" representation of SS'^y), y G is 

given by 

{6&(y), z)=e 2 Y,Yl d r V ( D %) ■ D rz(x), for z G 9. (48) 

xeif re.* 

Using bond integrals we rewrite this in a form that will be useful for our subsequent analysis. We 
will then use the bond density lemma whenever we need to transition between bond integrals 
and volume integrals. This process yields a notion of stress for atomistic models, which is 
related to the virial stress (see [2] for a recent reference; this connection will be discussed in 
detail elsewhere). A variant of this result for pair interactions in ID was developed in |32j . 

Proposition 11. Let y,z G 3 r , then 

(S^(y),z)= \T\Z a (y;T):Vz(T)= f £ a (y):Vzdx, 

Te.% n 

where the stress function £ a (y) G P*(^) 2x2 is defined as follows: 

E a (y;T):=^^- £ [d r V (D%y(x)) ® rU ' XT db. (49) 

Proof. For the sake of brevity we will write V x , r = d r V(D^gy(x)). 

Recall that T# = Ug62Z 2 (^ + T), for T G Sf e . It is easy to see that {xt* I ^ ^ ^e} is a 
partition of unity for M 2 . Hence, using the identity 

rx+er 

I V r zdb = D r z(x), (50) 

J X 



V r z db. 



we can rewrite (J48J) as 

Interchanging the order of summation and using the fact that V r z = (Vz)r holds db-a.e. (note 
that, if the bond is aligned with an element edge, then V r z is continuous across that edge) 
yields 

2 px-^sv 

(5^(y),z) =Y.\ T \Y.W\^ Vx > r 'T XT#Vz(T#)r db. 
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Finally, we use the fact that Dagy is 2Z -periodic to deduce that 
E a (y; T)= E 4r E ^ ® r ] f *T# 



db 



-2 



/' X~\~SV 

ref 1 1 zGJz? (,£2Z 2 x 

^2 rx+er 
Em ^ [ V ^r®r]l XT db. 

re,^ 5 1 1 xey# ^ 



(51) 



The term Vz(T # )r = Vz(T)r can be taken outside the bond integral, and hence, employing 
the identity 

a ■ (Gr) = (a® r) : G, for a, r G M d , G e M dxd , 

yields 

= E l T l E W\ [ V *,r®r]l XT# db \:Vz(T) 

--: E |T|S a (y;T): Vz(T). 

Te.% 



□ 



The terminology "stress function" for E a is motivated by the fact that £ a (y) takes precisely 
the same role as the first Piola-Kirchhoff stress tensor in the continuum theory of elasticity. 

In the next lemma we prove two useful properties of the atomistic stress function X a . We 
show that S a = dW under locally homogeneous deformations and give a quantitative estimate 
for the discrepancy between S a and dW . 



Lemma 12. The stress function E a defined in (49) satisfies 

£ a (y F ; T) = dW{F) VF e R 2x2 , T e 
Moreover, we have the estimate 

\^(y;T)-dW(Vy(T))\ < eM a osc(Vy; w|) Vj/GST, T G 17 £ , (52) 
where M a is defined in ^ and uj^ is defined in (42). 

Proof. Part 1: Since [d r V{F&) <g>r] is independent of x, we can apply the bond density lemma 
to the sum in curly brackets, to deduce that 

S a (y F ;T) = E [d r V{¥M) ® r|J JL £ I ' XT db = E [<^( 

Recalling that VF(F) = V(F&), it can be easily checked that the sum on the right-hand side of 
the second equality equals dW(F). 

Part 2: Let F = Vy(T), V x>r = d r V{D s y{x)), and V" F , r = d r V(F&). From part 1 we obtain 
that 

-2 r -] /»:r+er 



|£ a (y;T)-cW(F)| 



E XT 



db 



< 



Ei4 £ |i 



If 



XT db. 



(53) 
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We use the Lipschitz property [5] to estimate 

\V x , r - V F>r \ < Yl K,s\D s y(x) - Fa I 



X + ES 



V s y-Fs db 



< VM a Js| max IVy(x') - F|. 

— ' x' £(x,x+es) 

Since (x, x + es) C Wj., and recalling that F = Vy(T), we can further estimate 

max Vy(x') — F| < eosc(Vy; wf-). 

x'£(x,X+£s) 



(54) 



(55) 



We combine (|55j) with (54) and insert the resulting estimate into (|53j) to arrive at 

.2 rx+er 

I XT db. 



|£ a (y;T)-cW(F)| < eo sc(Vy; ^) ^ ^ |r| |s|M a s ^- ^ T 



An application of the bond density lemma, and referring to the definition of M a in ([6]), yields 
the stated result. □ 

6.2. The a/c stress function. We wish to derive a similar representation of 5(o ac in terms 
of a stress function S ac , as we did in £6.1 for bS & . A straightforward calculation along the same 
lines as the proof of Proposition 1 1 recalling first the definition of the partial derivative d^Ei 
from £2.3.6, yields the following result. 



Proposition 13. Suppose that ( |16| ) holds, then, for all y,z E <3f , 

(5^(y),z)= E \T\^ c (y;T):X7z(T)= j £ ac (y) : Vzdx, 

where 

S a (y;T), TG^ a , 
9W(Vy(T)), Te^, 



£ ac (y;T) := < 



E 



reM \T\ 



E, e ^# [^(^y(x)) ®r]f* +£r XT db 



^ +fr\E { x,x+er)^# [ 9 (x,x+er)Hy)®r]£ 
where \t ^ s a modified characteristic function: 



x+er ; 



1, x£dtt 



# 

i ' 



Xt(x), otherwise. 



Proof. Employing again the notation V x ^ r = d r V(D<%y(x)), the functional S£' SuC (y) can be written 



MS 



(5^(y),z)= j dW(Vy):Vzdx + e 2 J2J2 V ^' DrU ^ 

E d (x,x+er)Ei(y) ■ D r u(x). 
(x,x+er)£Si 



(56) 



The first term on the right-hand side of (56) gives rise to the definition of £ ac (y; T) for T G ^ c . 
(Note that (16) guarantees that the bonds in the second group in (56) do not contribute to O c .) 
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After the same calculation as in the proof of Proposition 1 1 , the second group on the right- 
hand side of (56) gives the definition of S ac (y; T) for T G ^ a , as well as the first group in the 
definition of £ ac (y; T) for T G 9*. 

The crucial modification to the previous argument is that the modified characteristic func- 
tions Xt> T ^ <3*e f° rm a partition of unity for J7i (except at a finite number of points, which do 
not contribute to bond integrals). Therefore, performing again a similar calculation as in the 
proof of Proposition 11 to "distribute" the third group on the right-hand side of (56) between 
interface elements only, we obtain the second group in the definition of £ ac (y;T) for T G ST*. 

Note that if we hadn't made the modification to the characteristic function, then S ac (y; T), 
T G ^ a U =^ c would contain contributions from E\. □ 



Since any discrete divergence-free tensor field may be added to S ac and still yield a valid 
representation of <5^ ac , it is not surprising that, in general, £ ac does not have the necessary 
property that £ ac (?/ F ;T) = dW{F) for all F G M 2x2 . This can already be observed in the nearest- 
neighbour, flat interface constructions in [44j . Hence, we need to construct an alternative stress 
function S ac representing e)<f ac that does have the desired properties. This construction will be 
undertaken in the remainder of this section, using the representation of discrete divergence-free 
vector fields as gradients of Crouzeix-Raviart functions discussed in §5.2| 



6.2.1. Consequences of global energy consistency. Recall from (13) that a functional $ G C {W) 
is called globally energy consistent if ^(yp) = ^(^f) for all matrices F G M 2x2 . The following 
lemma establishes a simple but crucial consequence of this property. 

Lemma 14. Suppose that $ G C 1 (^) is globally energy consistent, then 

{S*{W),VG) = N dW(F) : G V F, G G M 2x2 . 



Proof. From the assumption of global energy consistency, and (12), we obtain that 

<?(yF) = $M = \n\W(F) VFgM 2x2 . 
Since y? + tyc = yv+tG this implies that 

(S^),y G ) = |m lim W(F + tG)-W(F) = . Q □ 

N ' t^o t 

If we apply the foregoing lemma to an a/c functional <? ac we obtain the following corollary. 

Corollary 15. Suppose that <§ &c is globally energy consistent, then 

/ £ ac (y F )dx = cW(F) VF G M 2x2 . 
Jn 



Proof. This result follows simply from Lemma 14 and the fact that, for all G G 



d2x2 



(&&c(Z/f),Z/g) = / Sac(yF) : Gdx = ( / T, ac (y F ) dx ) : G. □ 
Jn \Jn J 
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6.2.2. Consequences of patch test consistency. The examples given in [44J show that patch test 
consistency does not necessarily imply that S ac (y F ;T) = dW(F). In the current paragraph we 
characterise the discrepancy between S ac (y F ;T) and dW(F). 



First, we show that the test functions Uh G ^ in the patch test (14) may be replaced by 
arbitrary displacements u£^, 

Lemma 16. Suppose that t§ &c is patch test consistent and that 2?^ U 5?? C 3? £ ; then we also 
have 

<^ac(y F ),u> =0 V^ef, FGl 2x2 . (57) 

Proof. Fix u£f ; then, using the assumption that ^ a U ^ C 2T e , we have 

Vudx = / u(g)fds= / If l u(&vds= I Vl^udx. 

n c J(dnf)nn J(dnf)nn Jn c 

Since S ac (y F ) = dW(F) is constant in O c , and since l^u = u in f2iUf2 a , we can therefore deduce 
that 

(£<£kc(yF),u) = / S ac (y F ) : Vudx = / S ac (y F ) : VI h udx = (S^dyp), I h u) = 0. 
Jo. Jn 

The penultimate equality requires some justification, but follows quite easily from the particular 



form of <p ac assumed in (18) and the assumption that 2?^ U^ 1 C 3T e . □ 



Lemma 17. Suppose that <§ &c is patch test consistent and globally energy consistent and that 
U STfi C 2T S ; then, for each F G M 2x2 , tfiere exists a function ip(F; •) G Nf (5J) 2 suc/i i/iai 

S ac (y F ;T) = aVF(F) + VV(F;T)J VTe^ 

where J is a rotation matrix defined in Lemma^ Moreover, j/int(fi a ) is connected, then we 
may choose tp(F) = in f2 a - 



Proof. If (f ac is patch test consistent then, according to Lemma 16 



/ S ac (y F ) : Vudx = {S^(y F ),u) =0 VnGt. 
Jf2 



Hence, according to Lemma r9l there exists a constant So G M 2x2 , a vector- valued Crouzeix- 
Raviart function tp = tp(F; •) G N^(j^) , and a rotation matrix J, such that 

s ac (y F ; t) = s + vv(T) J vr g sr £ . 

Using global energy consistency of <§ &c and Corollary [15] we obtain that 



8W(F) = I S ac (y F ) dx = S + / VVJ dx. 
Jn Jn 



If V-^dx = 0, then So = dW(F) and hence the result follows. 
To prove this, we integrate by parts separately in each element: 

/ VV>dx= ^2 ip®isds='^2 ® ^ + + ^~ ® v~) ds = 0, 

where, in the last equality, we used the fact that Jf{i> + — VO ds = for all edges /, since i\> is 
continuous in the edge midpoints. 



Finally, since S ac (y F ;T) = dW(F) for all T G ST^ (cf. Lemma 13 and Lemma 12), it follows 
that V0 = in O a . Hence, if int(f2 a ) is connected, then we can shift ip by a constant so that 
ip = in ft a . □ 
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6.2.3. The modified a/c stress function. We wish to construct a modified a/c stress function E ac 
that can be used to represent 5<f a c, and satisfies the crucial property that S ac (yp; T) = dW(F) 
for all F G M 2x2 , T G 

" to 



To this end, we generalize the Crouzeix-Raviart function VC~) defined in Lemma 17 



arbitrary deformations y G Since we will use later on that the modified function ip vanishes 
in O a , we require from now on that int(f2 a ) is connected so that we can choose V?(F, •) = in 
fi a for all F G M 2x2 . 

For each y G ^ and each face / G , / = T x n T 2 , we define the patch w/ = (Ti U T 2 ) \ ftf , 
and the deformation gradient averages 

F/(y) := 



£ f Vydx, if| W/ |>0 
0, otherwise. 



Note that was defined in such a way that ojf C cj^ n wt 2 . The value Fj(y) = for / C f2 a 
is of no importance, and could have been replaced by any other value. 

With this notation, and recalling the definitions of the edge midpoints qf and the periodic 

j± J I 

nodal basis functions from §|5.2.1| we can define 



kv\-)= E H F f(y)^f)cf- ( 58 ) 

Note, in particular, that ip(yf) = ^(F) for all F G M 2x2 . It is therefore natural to define the 
modified stress function 

S ac (y;T) :=£ ac (y;T)-W>(y;T)J, for T G (59) 

In the following lemma we establish some elementary properties of S ac . 

Lemma 18. Suppose that is energy and patch test consistent, that int(f2 a ) is connected, 
and that 2?^ U5JJ C 2? £ ; then the modified a/c stress function S ac , defined in (59), has the 
following properties: 

(5^ c (y),z)= / £ ac (y) : Vzdx Vy,zG^; (60) 
Jn 

S ac (y F ; T) = dW{V) V F G M 2x2 , ^; and (61) 

£ ac (y;T) = £ a (y;T) Vy G ^, Te^ a . (62) 



Proof. To prove (60) let z = ys + it for some B G M and «ef ; then 

E ac (y):Vz<Lc = (^ ac ( y ),3r)- / (V^J) : (B + Vu) dx. 



Since J^(VV'J) : Vu dx = by Lemma |9| and since f n S/tp dx = (see the proof of Lemma 17), 
the representation (60) follows. 

Property (61) follows from Lemma 17 and the fact that ip{y^) = ^(F): 

£ ac (y F ; T) = S ac (y F ; T) - V^(F; T)J = 3W(F) VT G 



Property (62) follows from the fact that we constructed tp{y) to be zero in fi a for all y G & , 
and from Proposition 13, which states that X ac (y;T) = £ a (y;T) for all T G J^ a . □ 
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6.3. The Lipschitz property. The final remaining ingredient for the proof of first-order 
consistency, is a Lipschitz property for £ ac , similar to the Lipschitz property (52) of E a . In 
order to ensure that there are no modelling error contributions from the atomistic region it 
turns out to be most convenient to work directly with the stress difference 



R(y;T) := S ac (y; T) - £ a (y; T). 
From ( 62 ) we immediately obtain that 

R(y;T) = VTe J/. 



(63) 



(64) 

In the remainder of the section we will estimate R(y;T) for T G 5* £ c U 3^. To motivate the 
following result we note that, from Lemma 12 and from (61) we see that 

1,2x2 



R(y F ;T) = VFg 



Tei c u X\ 



(65) 



hence, a suitable Lipschitz estimate for R yields the following result. 



Lemma 19. Suppose that all conditions of Lemma 18 hold and, in addition, that $\ satisfies 
the locality and scaling conditions (21) and (22); then 

\R(y;T)\<eM T osc(Vy;co T ) Vy e W, T £ 3%, (66) 



where Mt is defined in (45). 

The proof of this central lemma is split over the following paragraphs 
6.3.1. Estimates in the continuum region. Let T £ 3^ £ c , then 

| R(y; T) | = | [S ac (y; T) - V^(y; T)J] - £ a (y; T) 

< |dW(Vy)-£ a (y;T)| + |W(y;T)j| 

< eM a osc(Vy;uj T ) + |V^(y;T)|, 
where, in the last inequality, we used ([52]) and the fact that ojj, C wt for T E 3/~ £ . We still need 



(67) 



to estimate V?/>(y;T), which we postpone until £6.3.3 



6.3.2. Estimates in the interface region. Let T £ 3*^, and let V a 
R(y; T) = S ac (y; T) - S a (y; T) - V^(y; T)J 

1/ 



d r V(D M y(x)), then 



x<=JZ* 



® r 



f-x+er 2 r 



# 



bear 

b={x,x+er) 



rx+er 
J x 



xi-db 



EfflE K 



f Xr db-VV(y;T)J, 

J X 



which, after combining the first and third group, becomes 

rx+er 
J x 



6=(x,a;+er) 



2 t*X-\-&T 

Epfl E [^,r®^]/ XTdb- VVi(y;T)J 



=: R«(y;T)-R( 2 )(y;T)-V^(y;T)J. 
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We will again postpone the estimation of Vf/>(T) to f 6.3.3 and focus on the terms R^\y;T) 

and R (2) (y;^). 

Let F = Vy(T). Using the locality and scaling conditions (21) and (22), we can estimate 

J2 



f'X~\-ST 

RW( y] T)-RW{yr-T)\<^- £ \d b E { (y) - d^y^rU ' Xt 



db 



< 



fee^f 




(x,x+er) 




E 


E 






(x,x+er) ( x ' x 





M i r:S \D s y(x)-Fs\\r\J- XT db. (68) 



In the transition from the first to the second line we have used the fact that on bonds that lie 
on the boundary of fij the constant s is replaced by jM^, which effectively replaces Xt 
by xt- Bounding \D s y(x) — Fs| by the local oscillation, and applying the bond density lemma 
(see Lemma 12 for a similar calculation), we deduce that 

|R (1) (l/;T) - R (1) (y F ;T)| < £ MM-M*. osc(Vy; wr) = eMW(V 9 ; wt). (69) 

re.* se.^ 

Following closely the proof of ( |52[ ), we obtain a similar estimate for R^: 

|R( 2 )(y;T) - R (2) (y F ;T)| <£^^ \r\ \s\M* s osc(Vy; ur) = eM a osc(Vy; wr). (70) 

Note that it is enough to measure the oscillation over wp (which does not intersect with f2 a ), 
since (x, x + es) n fi a = for all s G \ -^ a # , se£ 

Combining (69) and (70), and using the fact that R(y F ;T) = 0, we conclude that 

|R(y;T)| = |R(y;T) - R(y F ;T)| 

< |R« (y;T) - R (1) (y F ;T)| + |R^(y;T) - R (2) (y F ;T)| + |VV>(y; T) — V^(y F ; T)| 

< e(M i + M a )osc(Vy; w T ) + | V^(y; T) - V^(F; T)\. (71) 



6.3.3. Estimates on ifi and on ip. To finalize the estimates in £6.3.1 and £6.3.2 we are left to 
establish a Lipschitz property for ip. The following result is a fundamental technical lemma 
that will allow us to achieve this. Its proof is deceptively simple, however, it uses implicitly 
many of the foregoing calculations. Moreover, some questions left open by Theorem 10 may be 
answered through a better understanding of this step. 

I I -ih. 

Lemma 20. Suppose that the conditions of Lemma 19 hold; then, for all f G , / C f2 c UOi, 
and for all F, G 6 M 2x2 , we have 

|V(F; q f ) - ip(G; q f )\ < e(M a + Afi)width(fii)|F - G|, 

where width(f2;) is defined in (46). 

Proof. Fix some /' G f C Q a ; then, for any connecting path 7 G ^f',f we have 
^(F; ff/ ) - V(G; 9/ ) = | (w(F) - W(G)) • dx 

= J (\dW(F) - E ac (y F )] - [dW(G) - E ac (y G )]) • dx. 
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Since E ac (y B ;T) = dW(B) for all Te J/U ST^ B G M 2x2 , the integrand vanishes in ft a U tt c . 
Hence, it follows that 

|V(F;g/)-^(G;g/)| < length( 7 n ^) max [£ ac (y F ; T) - £ ac (y G ; T)] - [0W(F) - dW(G)) 



Following closely the calculations in Section |6.3.2| we can deduce that 

hKF;g/)-V(G;g/)| < length( 7 n fii) (M a + M) |F-G|. 

Since we are free to choose the path 7, we can choose it so that length(7 n Oi) is minimized, 
which yields the stated result. □ 

Let T £ U ^ c , let F = Vy(T), and recall that are the Crouzeix-Raviart nodal basis 
functions associated with edge midpoints qy; then, using Lemma [20| we obtain 

|W>(y;T)-V^(F;T)| < £ |^(F/(i/); <?/) - V(F; g/)| | VC/(T)| 

/COT 



< (M a + Mi) width(Qi){ max |F/(y) - F|} je ^ |VC/(T)|1. 



/e 
/C9T 



/e*2 

/Cc>T 



A direct calculation yields 



£ |VC/(T)| = 2 + 2 + 2^^ 7. 



/CdT 

From the definitions of ff(y) and F it follows that 

|V^(y;T) - W(F;r)| < e7(M a + Afi)width(fii) osc(Vy; u T ) VT € 5" U 5; c . (72) 



Proof of Lemma 19. Combining (J72J) with (67) and noting that V?/>(F 
(61) and the fact that S ac (?/F, 



for T e 3T £ C (cf. 



arrive at the result of Lemma 



dW(F) in f2 c ), and also combining ( |72[ ) with (71), we finally 

□ 



6.4. Remarks on the conditions of Theorem 10 In this section, we construct simple ex 



amples to discuss the various assumptions of Theorem 10 We will show that most assumptions 
are also necessary. 



6.4.1. Technical conditions. The assumption that ^ a U Sfi C £F £ , and the assumption (16), 



were made for the sake of convenience of the analysis and simplicity of presentation. Dropping 
this assumption is not straightforward, but it is reasonable to expect that a careful analysis 
should allow to do so. 

The same statement applies to the assumptions made on the interaction potential; this was 
already discussed in §2.1.3 



6.4.2. Connectedness of£l a . The assumption that int(fi a ) is connected is more problematic; at 
this point it is unclear whether or not it can be removed in general. A more detailed analysis 
of the functions ^(F, •), F G M 2x2 , defined in £6.2.2 is required to understand this issue. There 



are, however, at least two special cases where one can attempt to remove it with relatively little 
effort: 

• Well separated components: If f2 a has several connected components, which are sepa- 
rated by an 0(1) distance, then one can localize the consistency error estimate to each 



of the components and obtain a qualitatively similar result as Theorem 10 
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71 







n 2 


72 



FIGURE 3. Atomistic region with two components to visualise the argument 



given in { 6.4.2 



Specific a/c methods: Suppose that the GCC method described in £2.3.4 is used to 
construct a patch test consistent coupling scheme, with parameters C X}TjS . Suppose, 
moreover, that Q a has two connected components, Q\ and O2, each of which have a 
portion of the boundary with the same orientation (say, normal ei), as displayed in 
Figure |3] 

It is then reasonable to assume that the parameters C X ^ )S have the same value in 
those parts of the interface surrounding and ^2 , which would imply that 

^ac(yF) ■ dx = S ac (y F ) • dx, 

'71 ^72 

Moreover, since S ac (?/p) — t^W^F) in the continuum region, we would, obtain that 



/ 



(£ ac (y F )-dVF(F)) • dx = 0, 



' 71U73U72 

where 7 2 denotes the curve 72 with reversed orientation. 

This shows that it is possible to choose -0(F; -) = in both components of f2 a , and as 



a consequence, Theorem 10 would remain true 



A related issue is the dependence of the modelling error estimate (44) on width(f2;), which 



comes solely from the Lipschitz estimate on F 4 ip(F;-); cf. £6.3.3 Hence, a better under 



standing of this function may also allow a finer analysis of this undesirable dependence. 
6.4.3. The global energy consistency condition. Global energy consistency is a natural and con- 



venient condition that yields the important intermediate result (Corollary 15) that 

Sac(yp) dx = dW(F) 



f 



VFG 



D 2x2 



(73) 



Note also that (73) implies ^(j/f) = ^(vf) + c for all F E 



02x2 



where c is a fixed constant 



that is independent of F; that is, (|73h is practically equivalent to global energy consistency. 



In some important situations patch test consistency already implies (73). The following 
result gives such a result for finite atomistic regions. 



Proposition 21. Suppose that S7 a U Oi C int(O); then patch test consistency (14) of 



implies (73). 



Proof. According to Lemmas 13 and |17| we have 

(<y<&c(i/F),u)= / dW(F):Vudx+ [ S ac (y F ) : Vudx. 
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Let G G M 2x2 and choose any u G ^ such that Vu = G in fii U $7 a ; this is possible due to the 
assumption that f2 a Uf2i C int(f2). Integrating by parts twice, letting v denote the unit outward 
normal to f2; U $7 a , and noting that the portions of the surface integrals along d£l cancel each 
other out, yields 



(<5<£kc(yF),«) 



I 8W(F) : Vudx+ [ £ ac (y F ):Gdx 

= -[ dW(F):(u®v)ds+[ £ ac (y F ):Gds 

ia(Q a uf2i) Jfi^ufii 

= -[ dW(F):Vudx+ [ S ac (y F ):Gds 

= / [S ac (y F ) - dW(F)] : Gdx. 

Since <o ac is patch test consistent, the last term vanishes, and hence the result follows. □ 

It turned out to be difficult to devise a counterexample, which clearly demonstrates that 
absence ( 73 ) can yield an inconsistent method. A more thorough investigation of this condition 
is still required. 



6.4.4. The locality condition. To show that the locality condition (21) (or a variant thereof) 
is necessary we assume, without loss of generality, that N is even and define a functional 



J(y) 



+ 



Y D eiV( X ] 



where, 



=5fi i± = {x G Jz? | xi < 0,x 2 = ±1/2}, 
Jz?j i± = {x G 5? | xi > 0, x 2 = ±1/2}. 



and 



From the definition it is obvious that J(y F ) = for all F G M 2x2 . Moreover, using summation 
by parts along the two lines JSfi + U + and JSfi _ U =5f}. _, it is easy to check that the patch 
test ( 14 ) holds. Finally, satisfies the scaling condition, 

21, if r = e\ and x, x' G Jz?^ b , a, 6 G {+,—}, 
0, otherwise. 



^(x,x+£r)^(x',a;'+£T-)^(y) 



However, ^ clearly violates the locality condition. 

We may think of & as an a/c functional for S & = 0, or, alternatively, as an additional 
contribution that can be added to any a/c functional whose interface satisfies Jz?i , U Jzfi _ U 

JSf| >+ U«5f|_ C int(fif). 

Let y G be "smooth" but not affine in the upper half plane {x G J2? | X2 > 0}, and let 
y = y& on the lower interface JSfi U then, testing 5 ^ {y) with the unique displacement 
-u G ^ such that 

u(x) = y(x) - Ax for x G {(0, 1/2), (1, 1/2)}; 

xi 1 — ^ u(x\, 1/2) is affine in [—1,0] and in [0, 1]; and 

D e2 u(x) = for all x G if ; 
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then we obtain, after a brief computation, 

Ha at Ml ^ (SJ(v),u) 



w: 



iVttl 



D ei y(x) - Aei 



x&se '+,+ 



Aei 



1/2 



This final estimate is scaled like a surface integral and is clearly a zeroth order term if y is 



smooth but y(Q, 1/2) 7^ y(l,l/2). This shows that the locality condition (21) (or a variant 
thereof) is indeed necessary to obtain a first-order consistency estimate. 

6.4.5. The scaling condition. It is fairly clear that the modelling error estimate can be arbitrar- 



ily large without the scaling condition (22). We nevertheless briefly discuss a simple example 
with a natural interpretation. 

Using a similar argument as in the previous paragraph, we define a functional ^ G C 2 (^), 
e 2 J, 

2 n \ -» In / \ 1 2 



J(y) = P 



\D ei y{x)\ \D ei y(x)\ , 



where /3 > is a constant, and where 

JS^ = {x € JSf |x 2 = ±1/2}. 

It is easy to see that J is patch test consistent, that J(yf) = for all F G M 2x2 , and that it 
satisfies the locality condition. 

Let y 6 *3f such that y(x) = Ax for X2 = —1/2, then testing 5 J? (y) with the unique 
displacement «£f such that 

u(x) = y(x) — Ax, for x G Jzf 1 , 

D e2 u(x) = 0, for all x G 

we obtain 



\^(y)\ 



w; 



> 



|| Vn|| L 2 



(3e 



e ^2 \D ei y(x) - Aei I 



nl/2 



(74) 



If y is "smooth" but not afiine, then the term in square brackets is of the order 0(1). By 
choosing /3 arbitrarily large, the modelling error can be made arbitrarily large as well. In 
particular, the choice j3 = 1/e would give a seemingly natural surface scaling to the interface 
functional, and in this case we would obtain an O(l) modelling error. 

7. Conclusion 

A fairly complete consistency analysis of general patch test consistent a/c coupling methods 
in (one and) two space dimensions was developed in this paper. The main result is the first 



order modelling error estimate, Theorem 10 The main undesirable condition is the assumption 
that int(fi a ) is connected. To remove this assumption a finer analysis of the corrector functions 
■0(F, •) defined in ^6.2.2 is required. At this point one cannot exclude the possibility that a/c 
methods exist for which this assumption is in fact necessary. 

Many open problems remain to be answered. First and foremost, one ought to answer the 



question whether a/c methods satisfying all the conditions of Theorem 10 always exist. In 
we present a general construction (a variant on the geometrically consistent coupling method 
|15j ) that appears to work in practise, however, we have no proof of this fact so far. Indeed, if 
it should turn out that in certain cases the "ghost forces" cannot be completely removed, then 



an extension of Theorem 10 estimating the contribution of the "ghost force" to the modelling 
error is highly desirable since such a result would provide the correct quantity that needs to be 
minimized. It is by no means clear that minimizing the "ghost force" itself is the best possible 
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target. A similar analysis would also be useful for estimating the modelling error of blending 
methods [55] , 

It should be conceptually straightforward, though technically more demanding, to generalize 
all results to higher order finite element methods in the continuum region, however, it would 
then also be desirable to obtain the second order modelling error estimate in the continuum 
region. Such a result seems difficult to obtain without a more detailed understanding of the 
corrector functions ^(F, •). 

An immediate question is whether a variant of the main result is still valid in 3D. This is by 
no means clear at this point. From a technical point of view, we require generalizations of the 
two main technical tools: the bond density lemma (£ |5.1[ ) and the characterisation of discrete 
divergence- free Po-tensor fields ( |5.2[ ). While the bond density lemma as stated in this paper 
is false in 3D, one can establish variants that are potentially usefull for a 3D analysis (work in 
progress). Generalising the explicit construction of {5.2 is entirely open at this point. 

Another important and difficult question is the extension to multi-lattices where the Cauchy- 
Born model is obtained through a homogenization procedure [TJ [5] . 

Finally, it should be stressed, that Theorem 10 is a general abstract result, and as such 
can undoubtedly be improved upon when a specific coupling method is analyzed. It may be 
possible for specific methods to obtain more information about the corrector functions tp(F, ■), 
and hence obtain a better estimate on the dependence of the modelling error on the interface 
width. For example, the consistency analysis in [II] requires no corrector functions at all, and 
in the consistency analysis of nearest-neighbour interactions [35] the corrector function vanishes 
in the continuum region. The proof of Theorem 10 may, however, serve as a general guidance 
for modelling error estimates in specific cases. 

Finally, the stability of a/c methods in 2D/3D is largely open at this point. 
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Appendix A. Proofs of ^3] 

Proof of Lemma^ For d = 1, since l e Vh — Vhi the result is trivial; hence assume that d = 2. 
Assume also that p < oo. Since the norms involved are effectively weighted F-norms, one can 
obtain the case p = oo as the limit p /* 00 • 

In this proof we will in fact use a periodic version of ( 40 ) , which is a simple consequence of 



i-x+er 



(gop (see also [UJ): 

\T\ =e 2 Y,i Xt* db. 

x<=J? Jx 

With our definition of the L p -norms for matrix-valued functions, we have 



2 2 



l v/ ^L ( n) = E / l v/ «l!> = E / l v ^ 

j=X j=X 



eyh\ p dx. 
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Using the periodic bond density lemma, and the fact that {xt# I T ^ ^e} is a partition of 
unity for M 2 , we have 



f \V ej I e y h \ p p dx = \T\\V e J E y h (T)\ p p 

Y Ki E y h (T)\yYf +£ej XT* 



db 



£ 



T&S? £ 
2 



x+eej 



\V ej I £ yh\ P p XT* db. 



|V ej / e y fe | p db. 



(75) 



We have also used the fact that V ej I e yh is continuous across edges that have direction ej. 

Due to the specific choice of the triangulation S? e it follows that V ej I e yh is constant along 
each bond (x, x + eej), and hence 



i 



x+eej 



|V^I £ y/i| p db = \D ej I £ y h (x)\ p 



L 



x+eej 



V e ,2/h db 



V rx+eej 
< 



px+eej 
J x 



where we employed Jensen's inequality in the last step. 



Inserting this estimate into ( 75 ) , and reversing the argument in ( 75 ) , we arrive at 

/•x+eej 



/ \V ej I e y h \ p p 6x < e 2 Y / \V ej y h \ P p db 
Jn xe^ Jx 



E^E/ K^|> # db 



E \t\hm t )\1 = \H. 



Vh\ 



□ 



Remark 9. From the foregoing proof, it follows that 

2 n 1/p 



Lp(O) 



□ 



Proof of Lemma^ To simplify the notation, we define the scalar function z = yi for some 
fixed i. Moreover, we prove the result only for d = 2; for d = 1 the result follows from the 
interpolation error estimates established in [42] . 

Step 1. W 2 '°° -interpolant: We first define a W 2, °°-interpolant z of z, using the C 1 - 
conforming Hsieh-Clough-Tocher (HCT) element |7j; see Figure [4j 

Let T G 2? £ and let Qt denote the set of vertices of T, and Ft the set of edges of T. 

For each vertex q G Qt-, we define the point value z(q) = z(q), and the gradient value by 
Vz{q) = j- ' c Vzdx, where 



\J{T' G {3r £ c )*\q£T'}. 



Similarly, for each edge / G Ft, f = T (IT' , with midpoint qy, we define the patch 
Slif n (T U T'), and the directional derivative Vvz(qf) = f^c V u zdx. 
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Figure 4. Illustration of the degrees of freedom in the C -conforming Clough- 
Tocher element: black dots denote point values, circles denote gradient values, 
arrows denote directional derivatives. 



Let <fq be the nodal basis function associated with the point value at a vertex q, <j>t the nodal 
basis function associated with the normal derivative at an edge /, and let &q t a be the nodal 
basis function associated with the a-component of the derivative Vz(g). 

Step 2. Estimating z — z: Fix T G x G T, and define F = Vz(T), then we have 

\Vz(T) - Vz(x)\ = |F- Vz(x)\ 



< 



feF T 



qeQr 
ae{l,2} 

Since all elements T G 3~ s are translated, scaled, and possibly reflected, copies of the reference 
triangle T = conv{(0,0), (1,0), (0, 1)}, it follows that the HCT nodal basis functions are given 
(up to translations and reflections) by 

tpf(x) = eiff(e~ 1 x), and & q , a {x) = e$ (?>a (e _1 x). 

Note in particular, that the gradients of these nodal basis functions are scale invariant, that is, 

||V99/||l°° < C and ||V$g )a ||L°° < C, 

where C is a fixed constant that is independent of e. 

From the construction of z it is easy to see that, for / G F q , q G Qt, ot G {1, 2}, 

: a ~ d Xa z(q) I < e osc( Vz; uj^); 



| Fuf — V Uf z(qj) I < e osc(Vz; w^), and F, 
and hence we obtain 

|| Vz(T) - V5|| LP(T) < Cie\T\ l ' p osc{Vz;uj^ 



(76) 



for some generic constant C\. 

Step 3. Interpolation error: Using standard interpolation error estimates [7], we obtain 



l^-VJ^|| LP(nc) <^||W 2 5|| LP(nc) . 



Let T G 3?s and F = Vz(T), then application of an inverse inequality, and ( |76[ ) yield 

|V 2 (z-z)|L m < C 2 e- 1 \\V~z-V\l p(T , <C 2 \T\ 1 /Po S c(Vz;co c T ). (77) 



IV 2 i| 



Lp(T) 



Finally, since z(x) = z(x) for all x G Jz?*, it follows that I^z = I^z, and hence we can 
estimate 



\Vz - VI h z\ 



LP (ft) 



\Vz - VI h z\ 



lp(o c ) 



< || Vz - V5-|| LP(nc) + ||Vz - Vhz\\ LP(nc y 



Employing (76) and (77), we obtain the stated result. 



□ 
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Proof of Lemma\4\ For each / £ let / = T_ nT+, T± £ let f-j- denote the corresponding 
unit outward normals, and uj'j = T + U T_ . 

We integrate by parts in each element ^ and use the fact that I £ Uh = u /i hi Vti U fi a to 
obtain 



^ / cr(T) : V(l £ u h - u h ) dx 
^ / ( a ( T +) u + + <t(T_)i/_) • (l E u h -u h ) ds 

< ^ eosc(a;uj' f ) / |j e «/, - u h | ds. 



Let v := /eU/j — u^. An application of [39, Lemma 6.6] yields the trace inequality 



(78) 



(79) 



IMb(/) < £ fIIlh^) + IIVuHli^). 

Furthermore, since v is Lipschitz continuous and v(p) = on every vertex of the triangulation 
^s, we can use [391 Lemma 6.8] to deduce that 

IMIli(w£) < V^e||V«|| L x( u ^. (80) 

Combining (80), (79), and (78), applying two Holder inequalities, and estimating the overlaps 
between the patches ui'j, we deduce that 



[&,I e U h ) - (<5>h,U h ) 



<(l + V2)e 08c(a)u/ f )\(J f \^ p \\Vi 



/ \ l/p 

< del \T\osc{a;oj^) p \ ||V 



Ilp'(Q c )- 



An application of Lemma [2] yields the stated result. 

Appendix B. List of Symbols 



□ 



a ■ b, a®b 
I ' I) I ' \p 

V r ,V 
dv 

V 



vector dot prod uct and tensor product; \ 
f-norms; qL3 



1.3 



weighted -P-norms; £1.3 



Lattice and lattice domain; \ 



2.1.1 



homogeneous deformation; £2.1.1 



spaces of periodic displacements and deformations \ 2.1.1 



periodic extension of a set or family of sets £2.1.1 



interaction range, £2.1.2 



finite difference operator and stencil £2.1.2 
directional derivative, deformation or displacement gradient, £1.3 



Jacobi matrix of vector valued function, £ 1.3 



atomistic energy, £2.1.2 



atomistic interaction potential, £2.1.2 
external potential in atomistic model, £j2.1.4 



5 2 <§', (■,■) first and second variations, abstract duality pairing, I 



d r V,d r ,V 



first and second partial derivatives of V, £2.1.3 



2.1.4 
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M r %,M a 
flT, h{x) 

Pi>Po>Pi ,P* 

ac 

^aj ^ci 
£y-a iy-c 

(x, x'), (x, x + er) 
Mi , M 1 



bounds on d r ^ s V and Lipschitz constant for 5^ a , £2.1.3 



triangulations, and edge sets, £2.3.1 



mesh size functions, i 2.3.1 



finite element spaces, £2.3.1 



finite element spaces, £2.3.1 



nodal interpolation operator for Pi(^) 
Cauchy-Born stored energy function, £ 2.3 



a/c energy and external potential, £2.3 



2.3.1 



atomistic, continuum, and interface region, £2.3.5 



atomistic, continuum, and interface triangulations, £2.3.5 



set of a tomist ic sites in a/c method, £2.3.5 



bonds, £2.3.5 



hi 


2.3.5 


[2.3.5 





scaled first partial derivatives of E- u £2.3.6 



bounds on second partial derivatives of E[, £ 2.3.6 



2?s i atomistic triangulation and edge sets, \ 3.1 



'S3 £ 

osc 



oscillation operator, £3.1 



4 



nodal interpolation operator for Pi(^), £3.1 



patch used in the interpolation error estimate, £3.2.2 



3.2.2 



Iw- 1 '* 
-i,p 



modified mesh size function { 
W^'-dual norm on Pi(5fc)*, £]3.4 

W 1,p '-dual norm on Pi(^)*, 93.5 



XniV n ,V n ,V n , V n 
XT 
XT 



<pi,4>2 first and second neighbour potential, §4 



notation for ID grid functions, ^4] 



characteristic function used in bond density lemma, { 5_A 
characteristic function used to define S ac , Prop. 



bond integrals, £5.1 
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P 



28 



Ni, Nf 

if, c f , cf 

f a ■ dx 

j 7 



Crouzeix-Raviart finite element spaces, £5.2.1 
midpoint of an edge / and associated nodal basis, £j5.2.1 
path integral 

Lemma |9[ p|23 



J rotation about 



TJT, 



width(fii) 

V x ,r, Vf,t 
Sac 

R(i/;T) 



atomistic interaction neighbourhoods, £j6| 
prefactors in modelling error estimate, Eq. (45), p 



width of Q h Eq. (46), p 



25 



25 



atomistic stress functionTEq. ( 49 ) , p 



26 



alternative notation for d r V(Dagy(x)) and for d r V(F&) 
a/c stress function, Prop, 
corrector function for E ac 
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2s 



P 

, Lemma 



17 



corrector function for S ac (y), Eq. (58), p 
modified a/c stress function, Eq. (59), p 
stress error, Eq. (63), p. 32 
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